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ABSTRACT 


This thesis addresses the problem of computing rapid slew maneuvers for 
a spacecraft antenna mounted on a double-axis gimbal with elastic joints. The 
performance of the system can be enhanced by designing antenna maneuvers in 
which the flexible effects are properly constrained, thus reducing the load on the 
spacecraft control system. The motion of a mass-spring-damper system is shown 
to be analogous to a spacecraft antenna slew with linear dynamics. This model is 
extended to a nonlinear double-gimbal mechanism with flexible joints, which 
better represents real spacecraft antenna dynamics. Rather than increase 
maneuver times to control flexible motion, this thesis presents optimal solutions 
that decrease maneuver times while allowing designers to easily constrain 
flexibility. Since it is impossible to recast the nonlinear system into a modal 
representation, an innovative approach is used to map the nonlinear dynamics 
into a linear system with a fictitious force. The fictitious force captures the effects 
of the nonlinearities so the vibrational motion can be constrained for a time- 
optimal slew. It is shown that by constructing an appropriate optimal control 
problem, the maneuver time for a flexible DGM can be decreased by 
approximately 42% compared to a conventional computed torque control 
solution. 
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I. INTRODUCTION 


NASA established the Tracking and Data Relay Satellite (TDRS) program 
in 1973 to provide continuous global services to critical missions in low-earth 
orbit while minimizing reliance on international ground stations. Since the first 
launch in 1983, NASA has commissioned and operated 11 TDRS in space (a 
twelfth satellite was destroyed in the Space Shuttle Challenger explosion in 
1986) [1], The current constellation consists of nine active TDRS at various 
positions in geosynchronous orbit (see Figure 1) to provide global relay capability 
for NASA missions such as the International Space Station and the Hubble 
Space Telescope [2], 
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Figure 1. TDRS constellation, from [2], 


Each TDRS uses multiple antennas to service its customers in an effort to 
reduce the number of satellites needed on orbit. Third-generation TDRS, for 
example, host a payload with two single access antennas and a multiple access 
antenna with 32 transmit and 15 receive elements [2], Each of these antennas 
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must be steerable in order to service moving targets in the low-earth orbital 
regime. This presents a unique problem in satellite design and operation in that 
the TDRS attitude control system must continuously accommodate the effects of 
slewing multiple antennas in order to maintain its orientation and minimize the 
associated disturbances on the spacecraft bus and other payloads. The 
actuators controlling the antenna gimbals must supply appropriate levels of 
torque to slew the antenna at a desired rate in addition to regulating flexible 
effects and mitigating pointing errors caused by the motion of other appendages. 
These functions are accomplished using various control modes. However, for 
slews, maneuvers are typically done slowly to avoid exciting the flexible modes of 
the satellite structures. In contrast, the spacecraft operators endeavor to 
maximize the utility of the TDRS by servicing as many customers as the control 
system will allow. Thus, the TDRS program faces a difficult problem in that 
demand for access is high and availability is limited due to the capabilities of the 
various satellite control systems [3], 

A. CONTROL METHODS FOR A FLEXIBLE STRUCTURE 

Several methods exist to enhance antenna control to increase user 
throughput. The antenna can be modeled as a manipulator arm connected to a 
controlling body through joints with some defined range of motion. The 
conventional method of computed torque control incorporates a simple approach 
for controlling the motion of such an arm [4], A proportional-derivative (PD) or 
proportional-integral-derivative (PID) feedback controller is common in computed 
torque control because these feedback laws can provide a near-linear response 
even for a system with nonlinear multi-body dynamics. A simple double integrator 
can then be used to model the motion of each joint of the arm for maneuver 
design. The feedback mechanism allows the controller to track the desired 
maneuver profiles [4], [5], This method of simplified maneuver design and 
implementation, while effective, typically dictates that the system be operated 
conservatively in order to reduce undesired vibrations in the antenna and other 
systems on a satellite. 
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Another method for antenna control applies Pontryagin’s principle to 
develop an optimal control based on the system dynamics and given boundary 
conditions [6], Pontryagin’s principle was applied in [3] and [5] to find minimum¬ 
time solutions to the TDRS antenna-slewing problem. The results show the 
potential for considerable savings in terms of time and increased customer 
access. 

For simplicity, most models, including the ones in [3], assume all bodies 
are rigid. Realistically, however, space structures typically exhibit some amount 
of flexibility. All materials possess physical properties that contribute to their 
response to forces, specifically impulses. The difficulty, however, lies in how to 
model this flexibility. Models typically require comparison to actual test results to 
ensure accuracy. As pointed out in [7], mathematical models help engineers 
understand the behavior of flexible systems in response to (or in the absence of) 
external input. These models, however, are usually simplified to avoid nonlinear 
dynamics because the nonlinearities cause unpredictable behavior, especially in 
the realm of flexible-body motion. Thus, modeling of flexible effects in mechanical 
systems is typically linear. This thesis explores the use of one such model to 
optimize antenna maneuvers while taking advantage of the physical properties of 
materials. 

B. THESIS OBJECTIVES AND SCOPE 

A significant challenge in designing control systems is to produce rapid 
slew maneuvers while minimizing the impact of flexible motion. An antenna must 
point at its targets with a certain degree of accuracy to ensure link closure and 
error-free transmission. A rapid slew is ineffective if it causes an antenna to 
“wobble” off target. Figure 2 provides a simple illustration of this phenomenon. 
The antenna or ground station must increase power to overcome the inaccuracy 
and maintain the link. To develop rapid slew maneuvers, the effects of flexibility 
must be considered. To this end, this thesis explores the concepts of modal 
analysis by looking at general linear mass-spring-damper systems with two and 


3 



three degrees of freedom (DOF). These linear systems are notional 
approximations of two-axis spacecraft antennas with linear dynamics. Applying 
modal analysis to these linear systems helps designers to understand and 
constrain vibrational effects of flexible systems. 



Figure 2. Illustration of a double-axis gimbaled antenna vibrating off 

target. 


Mapping modal analysis concepts to the antenna slew problem is difficult 
because the dynamics of the antenna motion are nonlinear. Hence, a goal of this 
thesis is to develop an approach for solving rapid antenna slews for nonlinear 
antennas that enables the flexible effects to be properly constrained in the modal 
space. 

Chapter II describes the basic concepts in modal analysis and discusses 
the value of frequency and impulse response functions to understanding the 
behavior of the flexible system. Two methods for producing a modal model are 
presented, one in the frequency domain and one in the Laplace domain. 
Chapter III expands the model to a 3-DOF system, representative of a notional 
spacecraft with an antenna mounted to a two-axis gimbal. Computed torque 
control is used to implement conventional slews to illustrate the behavior of 
flexible-body motion for the baseline system. In Chapter IV, optimal control 
theory is applied to minimize slew time in the context of several different 
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approaches for constraining flexible effects. Finally, in Chapter V, these results 
are applied to a more realistic spacecraft antenna model with nonlinear 
dynamics. Chapter V presents and illustrates an innovative approach for 
marrying modal analysis with nonlinear dynamical systems by treating the 
contributions of the nonlinear terms as a fictitious torque input. It is shown that 
the “new” system can be analyzed and optimized as a series of linear single-DOF 
systems but without the deleterious effects of a conventional linearization 
process. The thesis is brought to a close with some summary remarks and 
suggestions for future work in Chapter VI. 
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II. BACKGROUND ON MODAL ANALYSIS 


In modal analysis theory, the analysis of flexible effects is conducted in 
three phases. First, one must derive the governing equations of motion. That is, 
one must determine the values that populate the mass, damping, and stiffness 
matrices for the system. These second-order models are applicable because 
vibrational amplitudes are usually small and within the linear range of material 
properties. Second, one must perform a free vibration analysis of the system by 
setting any external forces to zero. The second phase of analysis produces a set 
of natural frequencies and damping ratios that lead to a matching set of mode 
shape vectors. These data are called the modal properties of the flexible 
system [7], The mode “shape” is defined as a unique “displacement pattern” [8] 
that describes the “relative displacements of all parts of the system” [7], Simply 
put, a vibration “mode” is one component of a set of ways by which a system 
oscillates in response to an input. The mode shape describes how each part 
moves relative to the others. Third, the forced response can be analyzed to 
understand the behavior of a flexible structure in response to an external 
input [7], which is of interest in many engineering systems. 

The general form of the governing equation of motion for a mass-spring- 
damper system with n degrees of freedom is 

Mx(r)+Cx(0 + Kx(0=f(0 (2.1) 

where MgM" is a diagonal matrix of mass values, CeM' !X " is a matrix of 
damping coefficients, KeR" M is a matrix of spring constants, f(?)el" !l is a 
vector of external forces acting on each mass, and x(r)eR' !Xl is a vector of the 
displacement of each DOF with respect to its equilibrium position. 


7 



Consider the two-mass system shown in Figure 3. 


Figure 3. 



c 
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*i(Q 

x 2 (t) 
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Example mass-spring-damper system with two degrees of 

freedom. 


The schematic describes a 2-DOF system in motion in a frictionless environment. 
The equations of motion for each mass are 


/j = m x x x + cx x + kx x - cx 2 - kx 2 


( 2 . 2 ) 


f 2 = m 2 x 2 - cx x - kx x + cx 2 + kx 2 . (2.3) 

Using (2.2) and (2.3), the mass, damping, and stiffness matrices in (2.1) may be 


assembled as M = 


m , 

0 


r -i 

i 


, c = 

c —c 

0 

m 2 


-c c 


, and K = 


k -k 
-k k 


The 


state and force vectors are x(r) = 


*i (0 

x 2 (t) 


and f(0 = 


m 

m 


. The 2-DOF 


system is used in this chapter to illustrate the different ways in which the modal 
properties of a flexible system can be evaluated. Undamped motion and motion 
with viscous damping will be studied in the absence and presence of external 
forces to understand the behavior of the system shown in Figure 3. 


A. NATURAL MODAL PROPERTIES 

A system without external excitation will exhibit unique modal properties. 

These properties are referred to as natural modal properties [7], The natural 

modal properties of a system describe the behavior of a forced system with 
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respect to certain frequencies. To find the natural modal properties, first consider 
the undamped case (c = 0) with f(?) = 0: 

Mx(f) + Kx(f) = 0. (2.4) 


Assume the solution to (2.4) is of the form x(t) = Xe imi , where XeM" xl 
denotes a vector of time-independent amplitudes. This leads to x(t) =-co 2 Xe im , 
which shows that the system oscillates following simple harmonic motion with 
frequencies co and amplitudes X[7], Substituting x(t) = Xe am into (2.4) gives the 
general eigensolution of the form 

-M(0 2 Xe iw + KXe iat =0 

(2.5) 

(K-A r M)X,=0 

where X r is the r-th eigenvalue. The corresponding eigenvectors, X r , are 
scaled versions of amplitudes X. For a system with n degrees of freedom, 
r = 1...« and X r =co 2 where co r is the r-th undamped natural frequency. 
Eigenvectors, v r , for each eigenvalue must be computed to determine the mode 
shapes of the system. The mode shape describes how the system responds at 
the corresponding natural frequency, and this shape is dominant at that 
frequency [9], In (2.5), (K-A r M)=0 must be true; otherwise, X r =0, the trivial 


solution. Let VeC' !Xl 
v = | 


be a matrix whose columns are eigenvectors such that 
I y ]. Thus, the modal model of a system is 


completely described by the two eigenmatrices 


and 



M- 


0 

0 

0 


0 

a 2 

0 

0 


o 0 

o 0 

0 



( 2 . 6 ) 


(2.7) 
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For the system given in Figure 3, it is desired to find the eigenvalues. 
Since this is a system with two degrees of freedom (n = 2) there are two modes, 
two eigenvalues, and two eigenvectors. Applying (2.5) gives 


k - X r m { 
-k 



k-X. 


m. 


X =0 


( 2 . 8 ) 


The eigenvalues, X r , may be solved as 


det 


k-X, 

-k 


m , 


k - Xm 2 


= m x m 2 X r - (m, + m 2 )kX r = 0 


(2.9) 


X x = = 0 

i „2 (m l + m 2 )k 
A>2 (Aj 2 


( 2 . 10 ) 


The eigenvectors can be found by substituting the eigenvalues back 
into (2.8) and solving for \ r (MATLAB’s eig command can do this symbolically). 

For A, = 0, 


k — / t,m, —k 

-k k-X 1 m 2 


1 

1 

s 

1 

0 


k -k 

0 


r * -i 

0 

—k k-X 1 m 2 

0 


-k k 

0 


o 

o 

0 


v„-v 21 =0 (2.11) 



v n 


v n 

a 

V 2 1 

= a 

v n 


where a is a scaling factor. Choosing a =— gives 


10 



( 2 . 12 ) 


For X n 


(m, +m 2 )k 


k - X 2 m x -k 

-k k- X 2 m 2 


v 2 =0 


k - X 2 m x -k 0 


-k k — X 2 m 2 0 


(2.13) 


, ( m,+m 2 )k 

k --— m. 


Equation (2.13) gives 


, (m, + m, )k 0 

k --— m 2 


(m { + m 2 )k 


(m x + m 2 )k 0 



i 2 

1 0 
m, 

1 0 

0 0 


m 2 

v 12 + —v 22 =0 
in, 


(2.14) 


\ 2 — cc 


a m. 
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Choosing a -—, 


22 


Therefore, 





ra, 


1 1 


(2.15) 


(2.16) 


In the first mode (described by the first column of V), both masses move in 
the same direction with the same amplitude, since both elements of v, have the 
same sign and value. The second mode is characterized by the masses moving 
in opposite directions due to the negative sign, and the relative amplitude of 

TYl 

motion is determined by the ratio of the masses (i.e., — [10]). Because the 

m t 

columns of (2.16) are eigenvectors, they are scalable in amplitude. However, the 
shapes of the modes, and hence the relative behavior of the two masses, will not 
change [7], The undamped free motion of the 2-DOF system is therefore 
completely described by 


[K] = [o 2 r ] = 


0 0 
0 (m l +m 2 )k 
m 1 m 2 


and 


V = 


1 


1 


’2h 

1 


(2.17) 


where \X r ] denotes a diagonal matrix of eigenvalues and [&r] denotes a 
diagonal matrix of natural frequencies, and r = l...n. 
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B. DIAGONALIZING THE MASS AND STIFFNESS MATRICES 


An important property called orthogonality can be used to diagonalize the 
mass and stiffness matrices in order to transform the original spatial model, 
which describes the motion in terms of displacements, x, to the so-called modal 
model, which describes the motion in a new set of coordinates. The 
transformation relies on the fact that 


M = V r MV 
K = V r KV 


(2.18) 


where Mel" is a diagonal matrix of modal masses, m r , and KgI“ is a 

_ 

diagonal matrix of modal stiffness values, k r , leading to A r = co; = — . The modal 

m r 

mass and modal stiffness matrices are not unique. Rather, they are based on 
how the eigenvectors are scaled. However, the ratio of these quantities is unique 
and fixed, giving unique modal frequencies, as in (2.17). The modal mass matrix 
is generally used to normalize the mode shape vectors. The normalized mode 
shape vectors can then be used to derive useful quantities such as effective 
mass or effective stiffness [7], The mass-normalized mode shape vector 0 f . eR' lxl 
is calculated as 


0 r = 



(2.19) 


such that column matrix OgR™ is 


<x> = 


0i 


02 



( 2 . 20 ) 


which has the property 

<&- 1 =(< b 1 ')- 1 ( 2 . 21 ) 

For the 2-DOF system given in Figure 3, the modal masses and stiffness 
values are 
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and 



1 

i 

r 

1 

, m 0 

_ 



m. 

0 

1 - 1 

M = y r MV = 

m 2 

i 



m. 


— 

0 

m 2 

1 




L 

- 1 

1 1 




1 

7 7U 

m l 

m 2 

L 




777, 

-m 2 

m 2 

1 

1 


m [ + m 2 0 




( 2 . 22 ) 




1 1 




K = 

V r KV = 

m 2 

k 

& 

1 - 1 

777 , 



- 1 1 

-k 

C 


1 



777 , 



1 1 



0 


0 1 

1 



= 

7 ( 7 ? 7 ,+ 777 2 ) 1 (m x +m 2 ) 

K /v 

1 

777 , 



777 , 

7 ? 7 | 

1 

1 



0 0 
0 

m x 


(2.23) 
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Using (2.22) and (2.23), it is possible to verify that [A r ] = M K as 
required: 


[A r ] = M K 


m x + m 0 

0 

-1 

0 

0 

0 

2 

m 2 

—- + m 2 


0 

k{m x + m 2 )” 


m 1 





m x + m 2 


m , 


m 2 (m, + m 2 ) 


0 0 
k[m x + m 2 )~ 


0 


m, 


0 0 
0 (Jjh + m 2 )k 
m 1 m 2 


(2.24) 


The mass-normalized mode shape vectors are 


o = 


01 02 


(m, + m 2 ) 


■i/ m. 


m, 


OT, 


N 1/ 
\/2 


^ (m, + m 2 )m 2 J 


(m, + m 2 ) /2 


in. 


X 1/ 

\/2 


(/»., + m 2 )m 2 J 


(2.25) 


Figure 4 and Figure 5 show an example of the relationship between the 
mode shapes and the displacement of each mass, for a sample system where 


m 1 =m 2 = k = l. The values of 


V2 V2 

2 2 


0.7071 0.7071 J (solid 


line) mean that each mass will move in phase and with the same amplitude at the 
frequency of the first mode, as seen in Figure 5(a). For the second mode, each 
mass will move in opposite directions but with the same magnitude. The motion 
is out of phase. In Figure 5(b), which shows the free motion of the 2-DOF 
system, mass 1 oscillates at the same amplitude as mass 2 but in the opposite 
direction. 
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Amplitude 


2 

1.5 

1 

0.5 

0 

- 0.5 

-1 

- 1.5 

-2 

Figure 4. Mass-normalized mode shapes for a 2-DOF system. 



Mass # 



0 2 4 6 8 10 

(a) 

Time (sec) 




Time (sec) 


Figure 5. Normalized displacements: (a) mode 1; (b) mode 2. 
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C. UNDAMPED CASE WITH AN EXTERNAL FORCE 

Adding a forcing function, f (t), to the 2-DOF system gives 

Mx(0+Kx(0=f(0- ( 2 . 26 ) 

The presence of the forcing function changes the procedure for modal analysis. 

In this case, the force vector can be written as f = Fe ,m where F e C" xl is a vector 
of complex amplitudes and co is the frequency [7], The variable co does not 
represent an output frequency as in (2.5); rather, it denotes an input frequency 
for the forcing function. The eigensolution in (2.5) then becomes 

(K-CQ 2 M)Xe iw =Fe im . (2.27) 

Following the procedure in [7], let 

X = H(£0)F (2.28) 


such that H e M' ,x " and 

H(ffl) = (K-ffl 2 Mr’. (2.29) 

In (2.29), matrix H is called a receptance frequency response function 
(FRF); it represents the system’s displacement in response to an input given in 
the frequency domain [7], Using the modal properties of the system, we may 
multiply by the matrix <t> and its transpose, as shown: 

O r H(«) '$ = O r (K-arM)® = O r K®-®WMO . (2.30) 

Since <t> is mass-normalized, the quantity O r KO gives a diagonal matrix of 

k 

-zr = (o; values. It also follows that ® r (y 2 M® = « 2 ® r M® = diag(or). 
m r 

Therefore, (2.30) may be rewritten as 

O r H(a») 1 0 = diag[(w 2 -« 2 )] (2.31) 

and further manipulated to give 

H(g))- 1 =(0 r r I diag[(« 2 -6) 2 )]0- 1 . (2.32) 
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Recalling that (AB) 1 = B 'A 1 , we obtain 


H(«) = 0(diag[(«;-(y 2 )]). (2.33) 

From (2.33), it is apparent that matrix H is symmetric. Because H 
represents the relationship between displacement and force, its components can 
be computed as 

H*(® ) = fr ( 2 - 34 ) 

F, 


where j = \...n and k = \...n. In other words, H x is the FRF for the 
displacement of the j'-th mass due to the k-\h element of the force vector, 
superposed over all the modes. The same relationship will be explored later in 
the Laplace domain. This “element-wise” relationship is useful to study the 
response due to certain forces on specific masses in the system [7], Equation 
(2.33) could also be used to find the full matrix of FRFs. 


H jk (co) 


I 


<E> O, 

jr kr 


co: 


CO 


" V V " R 

V Y jr y kr _ y r^jk 

ttm r (co;-co 2 ) ttco;-co 2 ' 


(2.35) 


The modal constant r R eC lxl is the specific receptance linking coordinates j 

and k for mode r [7], The notion of the modal constant, or “residue,” will also be 
seen later in the Laplace domain. 

Equation (2.35) highlights the principle of superposition with respect to 
modes—the total response of a system is the sum of the responses of each 
mode. In other words, 


H(g» = X,H«o) 

r =1 


where , H(co) denotes the FRF for mode r. In other words, 


(2.36) 


n n 

* II»► ■ 


r =1 k=\ 


(2.37) 
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For the 2-DOF system shown in Figure 3, the first mode dictates motion of 
both masses in the same direction and with the same amplitude. Thus, the 
receptance FRFs are the same. 


1 H 11 (©) = 


(m, + m 2 ) /2 (m l + m 2 ) 
-co 2 


-1 


or (m, + m 2 ) 


Using (2.38) gives 




,H n 

h 21 


!H 12 

1J 

l n 22 


-1 -1 

co 2 (m l + m 2 ) co 1 (m l + m 2 ) 
-1 -1 

co 2 (m l + m 2 ) ru 2 (m l + m 2 ) 


(2.38) 


(2.39) 


For the second mode, the amplitude of motion is the same for each mass, 
but the directions oppose each other. The symbolic representation for the second 
mode is cumbersome, so 2 H(ro) is not provided here. A numerical example is 

given later to clarify the details. 


D. VISCOUS DAMPING 

So far, only the undamped case has been studied. More realistic is the 
presence of damping in some form. Two general cases are discussed in [7]: 
structural (or hysteretic) damping and viscous damping. While realistic to some 
degree, there are advantages and disadvantages for each case. In this thesis, 
only viscous damping is considered. 

In the case of viscous damping, we first seek to determine the natural 
modes of the system by examining the unforced system, for which the general 
equation of motion is 

Mx(f) + Cx(f) + Kx(f) = 0. (2.40) 

The properties of each matrix are the same as previously defined for an n -DOF 
system. For the unforced case, with a solution of the form x = Xe h , where A = ico , 
the equation of motion yields 
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(MA 2 + CA + K)X = 0 . 


(2.41) 


The solution to this equation is a complex eigenproblem, with complex- 
conjugate-pairs of eigenvalues A and eigenvectors v such that 

(MA 2 +CA + K)v = 0 . (2.42) 

The bar symbol over the eigenvalue and eigenvector variables is used to 
distinguish them from the ones computed for the undamped free case. Each 

eigenvalue, A,., can be expressed as 

Vi 1 ??) (2-43) 

so that the eigensolution of (2.42) can be described as 

K’K 

and , (2.44) 

where r = l...n . In (2.44), the * notation is used to denote the complex conjugate, 
(i.e., X* r = a, (-^ r - ijl - C r ) )■ The orthogonality conditions are much more 
complicated than in the previous case. Let v be the eigenvector for the 
g-th mode; the same notation applies for A . Ewins [7] derives the following 
orthogonality conditions: 

(A r + A 9 )yjMv r + v fl Cy r = 0 

and (2.45) 

A,.A^v 2 Mv,. + v^Kv r = 0 
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If, and only if, v 9 = v* and X q = ft),. |-£ r - i-Jl - £, 2 j = A* then it is possible to 

determine the so-called natural frequencies co r and critical damping ratios of 
the system. 1 In this case, (2.45) leads to 


v H Cv 

2a)£ r = 

v, Mv. 


m.. 




and 

Kv, ^ 


v 

^7 » — 


Mv„ 


(2.46) 


where the superscript H denotes the Hermitian (complex conjugate) transpose. 
The similarity to the conditions in (2.18) is evident. 

The presence of damping makes the analysis of the forced system given 
in (2.47) more difficult to manipulate than for the undamped case. 

Mx(r)+Cx(r)+Kx(r) = f(r) . (2.47) 

To develop the FRF for the forced response, first recast (2.47) into state-space 
form [11], This can be accomplished using the following procedure if, and only if, 
the roots of the characteristic equation of the system have multiplicity one. Define 
state variable yeM 2 " xl as 


y = 


X 

X 


(2.48) 


so that the unforced equation of motion can be recast as 

I C M ]y + [ K 0 ]y = 0 . (2.49) 

This produces a system of n equations with 2 n unknowns, which is not helpful. 
Therefore, we make use of the identity equation 

| M 0 ]y + [ 0 -M ]y = 0 (2.50) 

1 This “natural frequency” matches the “undamped natural frequency” for single-DOF 
systems or systems exhibiting proportional viscous damping only [7]. 
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together with (2.49) to form a set of 2 n equations to give 


c 

M 

y + 

K 

0 

M 

0 

0 

-M 


This equation simplifies to canonical state-space form 

Ay + By = 0 


(2.51) 


(2.52) 


where AgR 2 " x2 " and BeR 2 ' 1 * 2 ' 1 are each square matrices. The eigensolution, 
assuming a solution to (2.52) of the form y = Y<? Xr , produces 2 n complex- 
conjugate-pair eigenvalues and 2 n complex-conjugate-pair eigenvectors 
(p r eC 2 ” xl . The eigenvalues and eigenvectors satisfy the general eigensolution 

(X,.A + B)(p r = 0 (2.53) 


where r = 1...2n. Orthogonality is useful here as well. Let ii) e C 2nx2n be a matrix 
whose columns are the eigenvectors, (p r , of the uncoupled system. Then, we 
may write 


\a r ] = $ T A& 
[b r \ = 


(2.54) 


such that 



(2.55) 


where a r and b r are the r-th elements of diagonal matrices [a,.]eC 2 ” 12 " 
and [z?,.] e C 2 " x2 ", respectively. It is important to note that the transpose operation 

in (2.54) is the non-conjugate transpose of the eigenvector matrix (versus the 
Hermitian transpose used in (2.46)) [7], 

The forcing function f can now be expressed as 


P = 


F 

0 


(2.56) 
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where PeC 2 " xl . To find the FRF, recall the procedure used in (2.28)-(2.35). The 
FRF can then be written as 


X 

icoX 


f <P?<P r _ 

r=i — A. r ) 


(2.57) 


which leads to the following individual receptance FRFs for the forced system 
(with viscous damping): 


H jk (co) 


s 


f 

VjAr 

, KK 


a r [a£ r +i((o-co r ^\-C t 

! )) a r (co£ r + i(cQ + 6) r jl-C: 

0 ), 


.(2.58) 


The relationship in (2.55) points to a useful coordinate transformation 
provided by the eigenvector matrix i3. Let y = $q where q e C 2 " xl is the so- 
called modal coordinate. From (2.53), one can see that q = e A ‘. As a result, (2.52) 
becomes 


A#q + B#q = P 


(2.59) 


The transformation by d puts the equations of motion in terms of modal 
coordinates to simplify further analysis. Premultiplying by yields 

tf r Atfq + tf r Btfq=tf r P . (2.60) 


So, for the r-th mode, (2.54) can be used to get 

^'A^q, + # r r B# r q r = t}'V = P' 
a r q + ^q,. = P; 


(2.61) 


Here, the r-th mode corresponds to the r-th element in the vector q [11], This 
transformation will be discussed further below. 
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E. PROPORTIONAL VISCOUS DAMPING 

A special case of viscous damping is called proportional (or classical) 
damping [12]. Proportional damping offers advantages in analysis over the more 
general case in that the mode shapes of a proportionally damped system are the 
same as those for the undamped system, and the natural frequencies are 
similar [7], The general equation of motion is the same as in (2.40), but now the 
damping matrix is proportional to the mass and stiffness matrices such that 

C = /3K + yM (2.62) 

where j8 and y are scalar constants. As a result, C is diagonalizable by the 
eigenvector matrix V to produce 

C = V 7 CV (2.63) 

where the r-th diagonal element of CgR'™ is the modal damping constant, c r , 
for the r-th mode. 


For the 2-DOF system in Figure 3, the modal damping matrix is 



1 1 

r 

1 

777 , 

c = v r cv = 


c -c 

-c c 

777, 


777, 


1 

1 


0 

0 

1 

777, 

(777, + 777, ) 

(777, + 777, ) 

777, 

-c- 

777, 

c- 

777, 

1 

1 


0 0 

(m, + m 2 ) 2 

U C 2 


(2.64) 


k 

because (3 = - exists and 7 = 0. 

c 

The modal transformation of a proportionally damped system is also more 
straightforward. The unforced system, from (2.18) and (2.63), is 
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(2.65) 


Mq + Cq + Kq = 0 

where q = V _1 x. Thus, the forced case is written in modal coordinates as 


Mq + Cq + Kq = V r f(0 = f'(0 ■ 


( 2 . 66 ) 


For the r-th mode, 


m A r + c r q r + k r q r =f r \t) 


(2.67) 


Equation (2.67) is clearly the equation of motion for a set of single-DOF mass¬ 
spring-damper systems, which represents the modes of a multiple-DOF system. 
Assuming a solution to the equation of the form q = Qe Af , where A is complex 
due to the presence of the damping terms, the r-th complex conjugate 
eigenvalue is computed from the characteristic polynomial 


m r X 2 r + c r \ + k r = 0 


A = ■ 


■ ±4 


—2 

C, 


■ 4m k„ 


2 m,. 




( 2 . 68 ) 


where, because of the proportional relationship in (2.62), 



and 





§G>r | 7 

2 2(0 r 


The damped natural frequency co' is written as 


or = co. 



We note that as long as is small, co' r ~ co r . 


(2.69) 


(2.70) 


For the 2-DOF system shown in Figure 3, by inspection, 
co 1 = (o[ = = Aj = 0, and for the second mode, 
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k(m ] + m 2 ) 


(2.71) 


co 2 = 


1 k 2 

k(m x + m 2 ) 7 

2 

m { 


2 

m 2 

—- + m 2 

m x 


C 2 = 


c(m, + m 2 y 

2 

in. 


c(m , + m 2 ) 


2 yfnA 


If 2 

\ 

f 

m. 




+ m 7 


' V m \ 

J 

V 


k(m ] + m 2 Y 

2 

m. 


\ 

2m ul 

r t 

< 2 \ 
ffl, 

—- + m 1 

) V 


l J 


(2.72) 


get 


For analysis of the forced response, it is possible to make use of (2.29) to 


H(ry) = (K + icoC - arM) ' 

V r H((y)V = V r (K + icoC - (y 2 M) ' V . 

Thus, the individual FRFs are represented by 


H *(®) = X 


V V 

jr kr 


k r + icoc r 


■co 2 m,. 


(2.73) 

(2.74) 


(2.75) 


where each term in the sum is the FRF corresponding to the r-th mode. 

For the 2-DOF system shown in Figure 3, the same symmetry exists in H 
when assuming proportional damping. For the first mode, r = 1, the FRFs are the 
same as for the undamped system (this may not always be the case) in (2.39). 

For r = 2, 


2 Hi,(ro) 


V V 

M2 y 12 


k 2 - icoc 2 - co 2 m 2 


( \ 2 

m 2 

V m x J 


k ( m ,+»!,) c ( m , + »!,) , 

V + 1(0 V ~( 0 ~ 


m: 


m. 


m: 


f my 

V m i 


+ m, 


k(m i + m 2 ) + icoc ( m , + m 2 ) - ft) (m ,ml + m x m 2 ) 
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(2.76) 



and so on. Numerical examples are discussed below that further illustrate the 
features of the modal analysis. 


F. NUMERICAL EXAMPLE 

Assume a nominal case where m l = m 2 =c=k = 1. For the undamped case 
with free vibration, 


w= 


0 

0 


0 

(m x + m 2 )k 
m,m 2 


0 0 
0 2 


(2.77) 


, m 2 



1-- 


1 -1 

m. 

= 


1 


1 1 

1 1 




It follows that ft),=0 and ca 2 = j2. Column vector v 2 indicates that the 
amplitudes of motion of both masses in mode 2 will be equal but in opposite 
directions. Thus, the modal matrices and mass-normalized mode shape vectors 
are 


with 



m l + m 2 
0 


0 



+ m 2 


2 0 
0 2 



k[m { + m 2 ) 


m. 


0 0 
0 4 


(2.78) 
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(2.79) 


o = 


(m, +m 2 )~ 
(m ] +m 2 ) 


m n 


m. 


m. 


N 1/ 
\/2 


(m, + m 2 )m 


m. 


2 y 
\^2 


(m, + m 2 ) m 


2 y 


V 2 V 2 
2 2 
V2 V2 
2 2 


Note that with different scaling of the eigenvectors in V, the modal matrices will 

contain different values, but the relationship A r = (£>)-— will be preserved. 

m r 

To analyze the frequency response in an undamped system with an 
external force, use (2.33) to get 


H(ffl) = 


V2 _ 

_V2 



-i 

V2 

V2 

2 

2 

-co 2 

0 


2 

2 

V2 

V2 

0 

2 -(o 2 




2 

2 




2 

2 

1 

1 

1 




2 co 2 

2 (or 

! - 2 ) 

ft ) 2 (ry 2 

- 

2) 



1 


1 


1 


(O 2 

(or-2) 

2 ft ) 2 2( 

ft) 

2 - 2) 



(2.80) 


Alternatively, each element of H can be computed individually using (2.35). For 
example, 


H l2 (co) 


'ji' 

'V2' 


' V2' 

's/2' 

2 J 

2 J 

_ 1 _ 

2 J 

2 J 


(2.81) 


A plot of the impulse response for the FRF in (2.81) is shown in 
Figure 6(a) as an oscillating linear growth toward infinity. Response to other 
inputs would show growth to infinity as well. This is not helpful in studying the 
oscillatory effects only, because the rigid-body motion of the system and the 
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oscillations are superposed, (i.e., the time history of Figure 6(a) is the sum of the 
rigid motion in Figure 6(b) and the oscillatory motion in Figure 6(c)). 



Time (sec) 

Figure 6. Undamped impulse responses in time domain: 
(a) H 12 — [H 12 + 2 H 12 , (b) ]H 12 , (c) 2 H 12 ■ 


For reference, the FRF for the first mode (from (2.81)) is 

I 

,H I2 (©) = —\ (2-82) 

CO 

and the FRF for the second mode is 
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1 

2 


(2.83) 


2 H 12 (ft>) - 

2-co 


Since the damping matrix of the system shown in Figure 3 is proportional 
to the stiffness and mass matrices (/3 = 1, 7 = 0), this system exhibits 
proportional damping and is straightforward to analyze. However, for 
completeness, the general complex modal case is computed to illustrate the 
steps. First, the unforced case is analyzed to compute the natural frequencies 
and damping ratios. First, we find the eigenvalues in (2.42). 


det(l 2 M +AC + k) = 0 


f 

- 

- 


- 

- 


r -A 


A 2 

0 


A 

-A 


1 -1 




+ 



+ 


V 

0 

A 2 


-A 

A 


-1 1 

L J y 


, A~ + A +1 —A — 1 
= det _ _ _ 

—A — 1 A’ + A+l 
= A r 4 + 2A r 3 + 2A r 2 = 0 
A, 4 + 2A, 3 + 2 A,. 2 = A r 2 (A,. 2 + 2A r + 2) = 0 


(2.84) 


Therefore, X i =X*=0. The quadratic formula yields A 2 =-l±/ from the second 
term. Following the same process in (2.11), for X i =X l *=0, we obtain the first 
eigenvector: 


A, + Aj +1 

-A, -1 

0 


1 

-1 

0 


1 -! 

0 

-A, -1 

A|~ + A| +1 

0 


-1 

1 

0 


o 

o 

0 


(2.85) 


For A 2 =-! + /, 
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(—1 + /) +(—l + z') + l 

-(-t + z)- 1 

0 

-(-t + z)- 1 

(—1 + z) +(—l + z') + l 

0 


—i -i 

0 


1 

1 

0 

-i -i 

0 


0 

0 

0 


v 2 = 

l 

±i 

0 

z 

-l 


0 


( 2 . 86 ) 


Recall that the eigenvectors occur in complex conjugate pairs, hence the ± sign. 

Use these eigenvectors and (2.46) to find the natural frequencies and 
critical damping ratios for each mode. For mode 1, 


®, 2 = 



_ v,“Cv, _ [ 

* i] 

i 

-1 

1 

Cl 

-i 

1 

L i 

zTz, 

Vj H Mv, 

[• *: 

i 

0 

1 



[ o 

1 

1 


r 


i 

-1 

1 

K 

_ v 1 w Kv 1 _ L 

1 1 ] 

-l 

1 

1 

zTz, 

v"Mv, 

[ 11 : 

i 

0 

1 



0 

1 

1 


= 0 


(2.87) 


which follows since A, = A* = 0. For mode 2, 


'2b 2 


= _£?_ 

_ \ 2 h C\ 2 _ ^ 

1 

->] 

1 

-1 

-1 

1 

1 

-1 

m 2 

v 2 "Mv 2 

[> 

-1 

1 

[ 0 

0 

1 

1 

-1 

_ k 2 

_ y 2 H Kv 2 _ ^ 

1 

-•]’ 

1 

-1 

-1 

1 

1 

-1 

m 2 

v 2 "Mv 2 

[> 

-1 

1 

0 

0 

1 

1 

-1 


= 1 = 2 
2 


= 2 


( 2 . 88 ) 


Therefore, a> 2 = V2 and £ 2 = . As a check, we may verify that these 

values reproduce the eigenvalues given by (2.84): 
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^ 2 =® 2 k +i >R 


A 2 = —1 —i 


2 [2 y 


: — 1 + i 


(2.89) 


Now we wish to derive the uncoupled modal motion of the system. The 
state-space form of this system (from (2.51)) is 


1 

-1 

1 

0 


1 

-1 

0 

0 

-1 

1 

0 

1 

y+ 

-1 

1 

0 

0 

1 

0 

0 

0 

0 

0 

-1 

0 

0 

1 

0 

0 


0 

0 

0 

-1 


y = o 


(2.90) 


Using the eigensolution form given in (2.53), the eigenvalues may be 
computed. The same eigenvalues as computed above should be obtained. The 
results are 


det[A,.A + B] = 0 


= det 


A r -A,. A r 0 

—A r A,. 0 A,. 

A, 0 0 0 

0 A. 0 0 


+ 


1-10 0 
-110 0 
0 0-10 
0 0 0 -1 


: det 


0 


(2.91) 


A ( .+ 1 -A, -1 A, 

-A,. -1 A,.+l 0 A,. 

A,. 0 -10 

0 A r 0 -1 

= A, 4 +2A,. 3 + 2A, 2 
= A,. 2 (X r 2 + 2A r + 2) = 0 

Again, A, = 0 and A 2 =-l±i. The eigenvectors, however, will be different 
because the system is now in the uncoupled state-space form. For A, = 0, 
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-l±i+l 

-(-1 + 0-1 

-l±/ 

0 

(-l±i)-l 

-1±/+1 

0 

-l±/ 

-1 ±i 

0 

-1 

0 

0 

-l±/ 

0 

-1 


±i 


-l±i 

0 

0 


±i 

0 

-l±/ 

0 

-i±/ 

0 

-1 

0 

0 

0 

-l±/ 

0 

-1 

0 


1 

0 

0 

1 1. 
— +—l 

0 

0 

1 

0 

2 2 

I±I,- 

0 

0 

0 

1 

2 2 

1 

0 

0 

0 

0 

0 

0 


1/ 


1/ 

72 


72 

-1/ 

72 

±i 

-1/ 

72 

-1 


0 

1 


0 


(2.93) 


Recall that the eigenvectors occur in complex-conjugate pairs, which is the case 


obtained in (2.93). The full matrix d is of the form 


(Pi Vl <Pl (Pi 


and produces matrix [a r ]eC 2nx2n of the form diag 


* * 


Ct i $2 Cly 0*2 


and 


[^ r ] — diag 


b { b 2 b* b* 2 


. From (2.54), 
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[a r ] = tf r Atf = 


r 

1 

1 

1 

-1 

1 

0 

1 

1 1. 

0 0 

-1 1 

1 _ 1 . 

-1 

1 

0 

1 

— ±—l 

-\- l 

1 

0 

0 

0 

2 2 

2 2 


0 

1 

0 

0 



-* 


0 0 

0 -2 


i i±i/ 

2 2 

i -I T I, 

2 2 


0 

0 


-1 

1 


(2.94) 




r 

1 


1 

-1 

0 

0 

1 

1 1 

0 0 

-1 1 

1 1 

-1 

1 

0 

0 

— ±— I 

— + —i 

0 

0 

-1 

0 

2 2 

2 2 


0 

0 

0 

-1 





0 

0 


0 

-2 + 2 i 


1 

1 

0 

0 


1 ±l,- 

2 2 

1 _ 1 
— + — i 
2 2 

-1 


(2.95) 


By inspection, it is evident that A 2 



= -l ±i, satisfying the relationship 


defined in (2.55). The FRFs for the forced system follow (2.58) for mode 2. Note 
that there is a division by zero in this expression for the first mode because of the 
double root at zero (multiplicity is two), which (2.58) cannot accommodate. 


It is evident that the computations for the general case of viscous damping 
are rather complicated. The computations become much simpler for a system 
exhibiting proportional damping. Because the example system studied here was 
specifically chosen to be proportionally damped, the eigenvectors are the same 
as for an undamped system, as shown previously. The damping matrix C can be 
transformed to the modal damping matrix of (2.64) if the following condition is 
true: 


(m- 1 k)(M" i C) = (m- 1 C)(m- i k) . 


(2.96) 


Evaluating (2.96) gives 
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(m _1 k)(m _1 c) 


(M“ 1 C)(M“‘k) 


1 

0 

-1 

1 

-1 

\ 

1 

0 

-1 

1 

-1 

\ 

2 

-2 

0 

1 


-1 

1 

J 

0 

V L 

1 


-1 

1 

/ 

-2 

2 

1 

0 

-1 

1 

-1 

\ 

C r 

1 

0 

-1 

1 

-1 

\ 

2 

-2 

0 

1 


-1 

1 

J 

0 

1 


-1 

1 

/ 

-2 

2 


.(2.97) 


Therefore, the modal damping matrix is 



C = 




0 0 
0 4 


So, for the unforced system, 


2 

0 

q+ 

0 

0 

q+ 

0 

0 

0 

2 

0 

4 

0 

4 


q = 0 


(2.98) 


(2.99) 


giving 


<7i = 0 

and 


q 2 + 2q 2 + 2q 2 0 . 


( 2 . 100 ) 


( 2 . 101 ) 


The eigenvalues are computed from the characteristic polynomial. For the first 
mode, by inspection (o l = C )l = X l = 0 (shown earlier). For the second mode, 



0 ), = fth 


Vl--CF = V2 1- 


V2 

v 2 y 


( 2 . 102 ) 


These values match those from the general case computed earlier. However, 
(2.102) shows that co 2 «.co 2 due to a high modal damping value relative to modal 
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stiffness. As discussed, the mode shape matrix V is the same as for the 
undamped system in (2.77). 

For comparison, the FRF for the first mass due to the second force is 
constructed below. 


H 12 (fl» 


V V V V 

_ ’ll ’21 _|_ ’12 ’22 _ 

k l + ia>c x - (0 2 m x k 2 + i(Oc 2 - (O 1 m 2 

_J_ 1 

-led 1 A + Aico-lo) 1 


(2.103) 


The impulse response to the first mode, by inspection, is linear towards infinity, 
so it provides no new insight. The presence of damping greatly changes the 
behavior of the second mode (compare Figure 7 with Figure 6(c)). The 
oscillations damp out quickly because the damping is large. Notice in Figure 7 
how the damping roughly halves the peak amplitude of oscillation of the 
undamped system. 



Figure 7. Proportionally damped impulse response of 2 H 12 for c = 1, 

c 2 = 4. 
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As an illustration of the effect of the relationship between the damping and 
stiffness constants, Figure 8 shows the impulse response of the second mode 
when the value of damping is decreased by an order of magnitude to c = 0.1, 
while keeping k = 1. Less damping increases the settling time of the system while 
also increasing its peak amplitude. 



Figure 8. Proportionally damped impulse response of 2 H 12 for c = 0.1, 

c 2 =0.4. 

G. FORCED VIBRATION WITH VISCOUS DAMPING IN THE LAPLACE 
DOMAIN 

Modal analysis can also be completed in the Laplace domain, where 
manipulation of the equations is easier. Many conventional control techniques 
have deep roots in the Laplace domain. To work in the Laplace domain, simply 
substitute co = into the FRFs discussed previously. 

Assuming zero initial conditions, Equation (2.1) describing the motion of 
the system in Figure 3 becomes Z(,y)X(,y) = F(s), where Z(s) = Mr +Cs + K. 

Matrix ZeM“" is called the dynamic stiffness matrix [9], The transfer function 
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HeM™ between the displacement XeM" xl and applied force FeM nxl is the 
inverse of the dynamic stiffness matrix and is given by 


H(s) = [Ms 2 + Cs + K]~‘ 


N(j) 
d(s ) 


adj[Mr+Cs + K] 
det[Ms 2 + Cs + K] ' 


(2.104) 


In (2.104), d(s ) is the characteristic polynomial, whose roots are the 

complex conjugate pole pairs, X r , of the system, and NeR"™ is a matrix of 
numerator polynomials that relates to residues (discussed below). The roots of 
the characteristic polynomial are the same as the eigenvalues of the chain of 
single-DOF systems, discussed previously. Recall from (2.68) that X r contains 
both a real part and imaginary part, such that X r = co£ r ± ico'. The roots of the 
numerator polynomials are the zeros of the system [13], 

The transfer function in (2.104) can be deconstructed into components 
separately describing motion of the system as a whole (termed “rigid-body”) and 
the oscillatory motion (“flexible-body”) using partial fraction expansion but only if 
the system, indeed, exhibits rigid-body motion. A more general case is discussed 
later. These components represent the modes from the previous section. For the 
2-DOF system shown in Figure 3, we have (from (2.104)) 


H(s) = 


m 2 s~ +cs + k 
cs + k 


cs + k 


m x s +cs + k 


s 2 (m t m 2 s° + (m, + m 2 ) cs + (m, + m 2 )kj 


(2.105) 


Using partial fraction expansion on each element of H(s) allows the 
individual transfer functions between the bodies to be determined. Let h,(i) 
equal the transfer function for the displacement X jk (s) and F*(s), such that 
X iJt (5') = H^(5)F t (5). In other words, X u (s) is the displacement of Xj due to /j 
and X 12 (s) is the displacement of x, due to f 2 . From (2.105), 
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H„(*) = - + 


5 2 m l m 2 s 2 + (m, + m 2 )cs + (/?/, +m 2 )k 

The numerator of (2.106) is 

a(m l m 2 s 2 +(m l +m 2 )cs + (m l + m 2 )fc)+ bs 2 = m 2 s 2 + C5 + /: 

Solving for a and ft we obtain 


(2.106) 


(2.107) 


a =- 

777, + 777, 

, 777.777, 

0 = 777,- 1 -— 

777, + 777, 


(2.108) 


Substituting (2.108) back into (2.106) gives 


777, 


H n 00 : 


■ + 


777, + 777, 


(777, + 777, ) S 2 777j777, 5 ? + (777, + 777, ) C.S' + (777j + 777, ) k 


(2.109) 


The first term in (2.109) represents the rigid-body motion of mass 1 (mode 1) 
while the second term represents the flexible-body motion (mode 2) due to its 
attachment to mass 2. Similarly, the remaining transfer functions are given by 


H,,(s) 


1 

( 777, + 777, ) .S' 2 


777,777, 

777, + 777, 

777,777,5” + (777, + 777, )C.S + (t 77, + 777, )/c 


H 21 (s) 


H,,(5) 


777,777, 

1 777, + 777, 

(?77| + 777, )5” 777,777,5 2 + (777, + 777, )CS + (777, + 777, )/c 


1 777, + 777, 

(777 1 + 777, ).V 2 777,777,5 2 + (777, + 777, )c5 + (777, + 777, )/t 


( 2 . 110 ) 


The symmetry is expected, as the motion of the two masses occurs only 
along the same axis, and they are only connected to each other. The first terms 
in each transfer function correspond to the rigid-body motion of the system. They 
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contribute to the first mode where co { =0. The second terms model flexible-body 
motion and represent the second mode. 

The residues of these functions are proportional to the modal vectors and 
can be used to construct the residue matrix R r eC“ [11]. (“Residue” is 
analogous to “modal constant” shown previously.) In general, assuming all poles 
have multiplicity of one [9], 


H(s) = £ 


r=l 


R 


R 


S-J 5 - A, 


V J 


( 2 . 111 ) 


where X r and X* are the complex conjugate pole pairs and n is the number of 
modes (or degrees of freedom). Residues are computed using partial fraction 
expansion, as above. The mode shape of mode r is described by a matrix 

'FeC“ whose columns are normalized modal vectors, y/ r , proportional to the 
corresponding elements of R,.. As was shown in the case of a proportionally 
damped system, the modal vector y/ r equals \ r from the undamped case. The 
modal vectors are calculated by using the residues from one row or column of H 
(for both modes). Interestingly, since the system is proportionally damped, the 
residues will be purely imaginary [11], 

The inverse Laplace transform of (2.111) gives the impulse response 
function h(t) studied earlier. 


hM=£(]*/■ 

r =1 



( 2 . 112 ) 


H. NUMERICAL EXAMPLE IN THE LAPLACE DOMAIN 


Assume a nominal case where m l =m 2 =c = k = 1. For the undamped 


case, A, = 0, A 2 = 2 , and V = 


-1 

1 


(from (2.10) and (2.16)), and it follows 


that the natural frequencies are co { = 0 and g> 2 = V2. Eigenvector v 2 indicates 
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that the amplitudes of motion of both masses in mode 2 will be equal but in 
opposite directions. 


The equations of motion in the Laplace domain become 

f 1 1 


X» = 
X 12 (s) = 
X 21 (s) = 
x 22 ( 5 ) = 


V5 

n 

2 


VS 


2 + r + 2.v + 2 


1 


s + 2s + 2 J 


Vf 

J 2 F ' 

'I F 

2 2 


1 


s + 2s + 2 


2 F| 


, 2 2 

V5 


1 


s + 2 s + 2 / 


-F, 


(2.113) 


It is readily apparent that the poles of the rigid-body component (r = l) are 
at s = 0. The rigid-body mode from the double-integrator terms in each equation 
implies that the system moves as a whole in the same direction due to the 
applied forces. The system’s second mode comes from the oscillatory motion 
between the two masses, modeled in the second terms. Take the transfer 
function component of the second term of X 12 (s). The poles and corresponding 

residues are 


X 2 = -\ + i 
2 R 12 = ± — 


(2.114) 


Recall that the poles match the eigenvalues calculated earlier. The residues can 
be used to reconstruct the transfer function from (2.111). 


H i2 (5) 


1 . 
—i 

4 


s -(-! + /) 


1 . 
—i 

4 


s-(-l-i) 


2 

s~ + 2s + 2 


(2.115) 


The residue matrix may be completed by following the same procedure. 
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(2.116) 


R, 


_ 1 . , 1 . 

+—i ±—i 

4 4 

, 1 • _ 1 . 

± — i + — i 

4 4 


The corresponding impulse response function is 


2 h l2 (t)-—ie^ ie ( 1 '^ = —— e 'sinf . 

2 12 4 4 2 


(2.117) 


The time history of (2.117) is precisely the same as the plot shown in Figure 7. 


I. MASS-SPRING-DAMPER SYSTEM FIXED TO A STRUCTURE 

Suppose now that the system from Figure 3 is connected on one side to a 
fixed structure, as in Figure 9. This connection represents a controller that may 
be implemented to stabilize the rigid motion of the flexible system. 



; m 


f 2 (t) 


= 0 


x 2 = 0 


Figure 9. Two-DOF mass-spring-damper system. 


The first mode of the system no longer exhibits rigid-body motion since the 
masses are fixed in place by the first spring-damper pair. The new equations of 
motion are 


/j = m,ic, + (c, + c 2 )i, + (&, +k 2 )x l -cx 2 -k 2 x 2 
f 2 = m 2 x 2 - c 2 x j - k 2 x x + c 2 x 2 + k 2 x 2 

which leads to 


(2.118) 
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. (2.119) 


Mx(f) + C x(t) + Kx(t) = f (t) 


in, 0 


C| + c 2 —c 0 


k t + k 2 — k 2 


o 

to 

x(0 + 


x(0 + 


x(r) = f(0 


1 

(N 

<N 

1 

_1 


—k 2 k 2 



Assuming the same nominal values, m l = m 2 = c, = c 2 = k { = k 2 = 1, the 
natural frequencies and mode shapes for the undamped, unforced system are 
given as 


M=[®r 2 ] = 


0.382 0 

0 2.618 


V = 


0.618 

1 


- 1.618 

1 


( 2 . 120 ) 


From (2.120), the modal matrices are obtained as 


0.618 1 

1 

0 


0.618 - 

- 1.618 



1.382 

0 

- 1.618 1 

_ 0 

1 


1 

1 



0 

3.618 

0.618 1 

2 

-1 


0.618 

- 1.618 



0.528 

0 

- 1.618 1 

-1 

1 


1 

1 



0 

9.472 

0.618 1 

2 

-1 


0.618 

- 1.618 



0.528 

0 

- 1.618 1 

-1 

1 


1 

1 



0 

9.472 


( 2 . 121 ) 


The equation of motion in modal coordinates is, therefore, 


1.382 

0 

q+ 

0.528 

0 

q+ 

0.528 

0 

q = 

0.618 1 

0 

3.618 

0 

9.472 

0 

9.472 


1 

t—H 

00 

t—H 

VO 

r-H 

1 

_ 1 


The equations for each mode, after dividing through by m r , are thus 

< 7 , + 0 . 382 ^ + 0 . 382 ?1 = 0 . 447 /, + 0 . 724/ 2 = (2.122) 

1.3o Z 

q 2 + 2 . 6184 , + 2 . 618 ^ =- 0 . 447 /, + 0 . 276 /, = -^~ . (2.123) 

3.618 


These equations describe the motion of two single-DOF systems, where the 
motions of the modal masses are uncoupled and are of the form 
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q r + 2C r CQ,.q r + CO r 2 q, 


(2.124) 



For the forced system, 


co l = V 0.382 = 0.618 
0.382 


c 


(2)(0.618) 


0.309 


©,' = 0.588 


and 


(0 2 

c 2 


= a/2.618 = 1.618 
2.618 


(2)(1.618) 


= 0.809 


ffl'=0.951 


(2.125) 


(2.126) 


Equations (2.125) and (2.126) can be interpreted as follows: mode 1, with 
the smaller natural frequency, represents the main motion of the system, 
whereas mode 2, with the larger natural frequency, represents the flexible effect. 
This flexible effect perturbs the main motion of the system. In this context, it is 
desirable to minimize the magnitude of q 2 and its derivatives since the flexibility 
is parasitic to the performance of the control system. This result illustrates the 
primary reason why modal analysis is a valuable tool in control system design. 

The controlled 2-DOF system shown in Figure 9 cannot be deconstructed 
into rigid-body and flexible-body components, as shown in the previous section, 
since motion of both masses is oscillatory. Assuming the same nominal values 
m x = m 2 = c, = c 2 = k x = k 2 = 1, in the Laplace domain the FRF is 



S 2 +5 + 1 

5 + 1 

N(s) adj[M,r +CX + K.] 

5 + 1 

5’ + 2 5 + 2 

~ d(s) ~ det[Mr +Cs + K] 

5 4 + 35 3 +45 2 +25 + l 


(2.127) 


where the poles are X, =-0.191 ±0.588/ and X 2 = -1.311 ±0.951/', which match 
the values computed previously. The residue matrices for each mode are 
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(2.128) 


+ 0 . 235 / + 0 . 380 / 
+ 0 . 380 / + 0 . 615 / 

+ 0 . 380 / ± 0 . 235 / 
± 0 . 235 / + 0 . 145 / 


Therefore, the mode shape matrix, using the first row of each residue matrix, is 


R 11 

2 r„ 

r 2 , 

2^21 

r 2 , 

D 

2^21 

r 2 , 

D 

2^21 


- 0.235 - 0.380 

- 0.380 0.235 

1 1 


0.618 - 1.618 

1 1 


(2.129) 


While determining mode shapes and residues in the Laplace domain is 
straightforward, the transfer functions in H cannot be used to get the modal 
coordinates of the system. They do, however, provide meaningful 
representations of the behavior of the system in physical space. 


J. SUMMARY 

The motion of flexible multi-body systems is complicated due to the 
interaction between the various degrees of freedom. Modal analysis enables the 
description of the motion to be decomposed into a series of single-DOF systems, 
each with different modal properties (i.e., frequency and damping). Control 
applications are typically concerned with manipulating the behavior of the 
dominant (slow) modes while minimizing the excitation of the flexible (fast) 
modes. Thus, modal coordinates are a natural space for synthesis of rapid slew 
maneuvers, and this chapter has provided the necessary background for 
accomplishing this task. 
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III. CONVENTIONAL CONTROL OF A 3-DOF SYSTEM 


This chapter increases the complexity of the mass-spring-damper system 
by adding a third DOF. The system is now analogous to a controlled spacecraft 
body with a flexible antenna mounted on a double-axis gimbal. The controller is 
used to maintain the spacecraft’s position and move the links of the notional 
antenna. This chapter analyzes a model to understand the behavior of such a 
system from the perspective of conventional control techniques. 

A. DEFINING THE MODEL 

Consider the three-mass system shown in Figure 10. It is representative of 
a controlled spacecraft with a 2-DOF flexible appendage. In this case, the 
stiffness k { and damping c, can be attributed to the feedback control of the base 
body, whereas k 2 , k 3 , c 2 , and c 3 are the stiffness and damping constants of the 
flexible appendage (locked in a specific orientation). 



x^ = o x 2 - o x 3 = o 


Figure 10. Mass-spring-damper system with three degrees of freedom. 

The equations of motion governing this system are 

/ = m 3 x j + (q + c 2 ) i , — c 2 x 2 + (k x + k 2 )x , - k 2 x 2 

f 2 = m 2 x 2 - c 2 x j +(c 2 + c 3 )x 2 - c 3 x 3 - k 2 x ] + ( k 2 + k 3 )x 2 - k 3 x 3 . (3.1) 

f 3 = m 3 x 3 - c 3 x 2 + c 3 x 3 - k 3 x 2 + k 3 x 3 
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Thus, 


m l 

0 

0 


c,+c 2 

-c 2 

0 

0 

m 2 

0 

C = 

-c 2 

C 2 +C 3 

-c 3 

0 

0 

m 3 


0 

-c 3 

c 3 


k j + kj 

—k 2 

0 


*i (0 


' f(0 ’ 


-k 2 

k 2 + k 3 

-K 

- x(/) = 

x 2 (t) 

, and f (0 = 

/ 2 (0 

. Unit values for 

0 

~ ^3 



x 3 (t) 


I 

i_ 



each parameter in the system are assumed for convenience so that 

m 1 = m 2 = m 3 = c, = c 2 = c 3 = k { = k 2 = k 3 = 1. 


B. ANALYSIS OF THE COUPLED SYSTEM 

Using the method described in Chapter II, the transfer function matrix for 
the system in Figure 10 can be written as 


H (s) 


s 4 + 3 s~ + 4 s~ + 2 s + 1 5' + 2 s~ + 2s + 1 s~ + 2s + 1 

s 3 + 2s 2 +2s + 1 s 4 + 3s 3 + 5s 2 + 4s + 2 s 3 + 3s 2 + 4s + 2 
s 2 + 2s +1 s 3 + 3s 2 + 4s + 2 s 4 +4s 3 + 7s 2 + 6s + 3 

s 6 + 5s 5 +1 Is 4 + 13s 3 + 9s 2 + 3s +1 


(3.2) 


Recall that the poles of this transfer function matrix are the roots of the 
denominator polynomial (or “characteristic equation”), which are equal to the 
eigenvalues of the system 

[A,] = diag[ -0.990 + 0.4339/ -0.7775 ±0.9749/ -1.6235 ±0.7818/ ]. (3.3) 


Using the procedure discussed in Chapter II, the undamped mode shape vectors, 
rounded to four decimal places, are 



1 1 1 
1.8019 0.4450 -1.2470 

2.2470 -0.8019 0.5550 


(3.4) 


The modal matrices, then, are 
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(3.5) 


M = V r MV = diag[ 9.2959 1.8412 2.8629 | 

C = V r CV = diag[ 1.8412 2.8629 9.2959 ] . 

K = V r KV = diag[ 1.8412 2.8629 9.2959 j 

From these equations, the natural and damped frequencies, co r and co' (in 
radians per second), as well as the critical damping ratios, , may be computed 
as 

[«,.] = diag[ 0.4450 1.2470 1.8019 j 

[£,.] = diag[ 0.2225 0.9749 0.7818 ] . (3.6) 

[fi>;] = diag[ 0.4339 0.9749 0.7818 j 

The relative values of the modal masses in (3.5) help to illustrate to an 
important property of this system. Modal mass m 1 = 9.2959 is five times larger 
than rh 2 and over three times larger than m 3 . The damping ratios increase 
through the modes. One could hypothesize, then, that the motion of mass 1 
dominates the total motion of the system, just as was observed in the system 
with rigid-body motion (see Chapter II). Analysis of the FRFs and IRFs will 
support this hypothesis. Figure 11 shows the Bode frequency response and 

Figure 12 shows the motion of mass 1 in response to impulse force 1 for all 

modes. 
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Figure 11. Bode magnitude and phase plot of H n (to) 

for a 3-DOF system. 


In Figure 11, the vertical dashed and dot-dashed lines represent the damped and 
natural frequencies for each mode, respectively. Clearly, the total motion of mass 
1 due to force 1 is greatly impacted by mode 1, especially in the low frequency 
range. This is because the magnitudes of mode 2 and mode 3 are 
-10 dB and -20 dB, respectively, in the low frequency range. In the high 
frequency range, the responses of modes 2 and 3 dominate. This is why flexible 
systems oscillate; sudden changes in the applied force excite these high 
frequency effects. The impulse responses (Figure 12) confirm that mode 1 
dominates the response and that modes 2 and 3 decay quickly due to the large 
damping. If this damping is reduced, the higher modes will be superposed onto 
mode 1 to a greater degree. This is the characteristic that degrades the 
performance of a control system and should be avoided via a combination of 
control and reference trajectory design. Figure 13 illustrates this point; it shows 
the IRF h n (r) for reduced damping (c, = c 2 = c 3 = O.l). 
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Figure 12. Impulse response function h n (r) for a 3-DOF system with 

“heavy” damping (cj = c 2 = c 3 = l). 



Figure 13. Impulse response function h u (r) for a 3-DOF system with 

“light” damping (c, = c 2 = c 3 = O.l). 
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Figure 14 shows the impulse responses for the transfer function matrix in 
(3.2). The repeated impulse responses (as for H 12 and H 2I , for example) result 
from the symmetry of the matrix H(s). It is apparent from the figure that m 2 and 
m 3 experience the greatest displacement about their equilibrium positions due to 
the forces acting on them. This follows due to the nature of the system; each 
mass is less constrained than its previous counterpart as one moves away from 
the fixed structure (base). In the control of a multi-body system, this is important 
because if forces 2 and 3 are not carefully controlled, they can cause large 
displacement of the base body (i.e., mass 1). See responses for H l2 and H l3 in 
Figure 14. 



Figure 14. Impulse responses, h jk (t), for system in Figure 10. 
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The residue matrices for each mode, using partial fraction expansion from 
Chapter II, are 

+0.1240/ +0.2234/ +0.2786/ 

R, = +0.2234/ +0.4025/ +0.5019/ 

+0.2786/ +0.5019/ +0.6259/ 

+0.2786/ +0.1240/ ±0.2234/ 

R, = +0.1240/ +0.0552/ ±0.0994/ | . (3.7) 

±0.2234/ ±0.0994/ +0.1791/ 

+0.2234/ ±0.2786/ +0.1240/ 

R 3 = ±0.2786/ +0.3473/ ±0.1546/ 

+0.1240/ ±0.1546/ +0.0688/ 

These residues may be used to construct mode shape vectors 

D D D 

1 2 XV 11 3 XV 11 

1 R11 2^11 3^11 

D D D 

Vp _ 1 xx 21 2 r *'21 3 rv 21 

1 R11 2^11 3^11 

|i n i> 

1 xx 31 2 .31 3 rv 31 

D D D 

1 AV 11 2 iv ll 3 ix -ll 

1 1 1 
= 1.8019 0.4450 -1.2470 

2.2470 -0.8019 0.5550 

Figure 15 plots these mode shapes. Recall that a plot like Figure 15 is useful in 
understanding the relative motions of each degree of freedom of the system for 
each modal frequency. Each line in Figure 15 represents a single mode mapped 
to each DOF. For instance, motion of mode 1 (blue line) causes all three masses 
to move in the same direction (because each amplitude has the same sign) but 
with varying magnitudes. One should expect m 3 to oscillate at more than double 
the amplitude of m l (2.25:1 ratio) and only slightly more than m 2 (1.80:1 ratio). 
Mode 2 (red line) would see m 3 moving in a direction opposite the other two 


1 1 1 

-0.2234/ -0.1240/ 0.2786/ 

-0.1240/ -0.2786/ -0.2234/ 

-0.2786/ 0.2234/ -0.1240/ 

-0.1240/ -0.2786/ -0.2234/ 
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masses because its amplitude has the opposite sign of the other two, and so on. 
(It can also be shown that these mode shapes can be formed from the peak 
amplitudes of the imaginary parts of their FRFs.) 



Figure 15. Mode shapes of 3-DOF system (c, = c 2 = c 3 = l). 

C. MODELING A SPACECRAFT SLEW 


With minor modifications, Figure 10 can be used as a representative 
model for a spacecraft antenna slew. Suppose a steerable antenna is attached to 
the body of a spacecraft (see Figure 16) by a two-axis gimbal. For simplicity, 
assume the center of rotation of the gimbal is aligned with its center of mass, and 
the antenna is mounted in such a way that there is no offset between it and the 
spacecraft body. 
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Figure 16. Notional schematic for a two-axis antenna mounted on a 

spacecraft body, after [5], 

The system in Figure 10 in its current form can be used as a notional 
model for the flexible component of motion (vibrations) resulting from the forces 
applied to the azimuth and elevation bodies, where the spacecraft body is 
analogous to m x , azimuth to m 2 , and elevation to m 3 . The spring and damper 
connecting mass 1 to the rigid structure represents the spacecraft controller that 
seeks to counteract the motion induced by the antenna and maintain the 
spacecraft in a certain position. The forces applied to the antenna to complete 
the slew are the forces that are input into the equations of motion in (3.1). The 
forces are also the inputs to simulations of the transfer function H(s) in (3.2). 

In order to model the entire maneuver—not only the flexible motion in the 
notional antenna’s locked position—Figure 10 must be modified to include the 
gross motion of the notional gimbals. A revised model is shown in Figure 17. The 
model shows that masses 2 and 3 have moved to new equilibrium positions 
defined by the desired final position at time t f , while the flexible-body motion has 

been preserved and referenced to a non-stationary equilibrium point. 
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*l(ff) ® X 2 (ff) - X 2 f *3(ff) *3f 

(b) 

Figure 17. Three-DOF system modeling rigid translation with flexible 


effects: (a) initial position; (b) final position. 

Let X 2 = x 2 (t)-x 2 (t Q ) and X 3 = x 3 (t)-x 3 (t 0 ) be the translation of each mass as the 
prismatic joints are moved. The motion of the masses now also includes what 
can be equated to rigid-body motion. The total motion of the masses is the sum 
of this rigid-body motion and the flexible-body component for each mass. Thus, 
the equations for the total motion become 

/, = m x X x + m l x i + (q + c 2 )xj - c 2 x 2 + (k ] +k 2 )x 1 - k 2 x 2 

f 2 - m 2 X 2 + m 2 x 2 - c 2 i, + (c 2 + r 3 )i 2 - c 3 i 3 - k 2 x x + (k 2 + k 3 )x 2 - k 3 x 3 (3.9) 
f 3 — m 3 X 3 + m 3 x 3 — c 3 x 2 + c 3 x 3 — k 3 x 2 + k 3 x 3 

where X x =X 1 = X 1 =0 since the translation of the base is zero. The total motion 
of the system is the linear superposition of the rigid-body motion and the 
vibrational effects. 
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D. 


SLEWING WITH COMPUTED TORQUE CONTROL 


Assume the following boundary conditions and constraints for a notional 
antenna slew are: 

(x 1 (r 0 ),X 2 (r 0 ),X3(r 0 ),X 1 (r 0 ),X 2 (r 0 ),X 3 (r 0 )) = (o,-lO,5,0,0,0.l) 
(x 1 (r / ),x 2 (r / ),x 3 (r / ),x 1 (r / ),x 2 (r / ),x 3 (r / )) = (0,l5,-l5,0,0.l,0) 

|* m J = 0.75 , (3.10) 

|/J = l-5 
[t 0 ,t f ) = (0,38.1) 

where X 2 and X 2 are the translation and velocity, respectively, in the notional 

azimuth plane, and X 3 and X 2 are the translation and velocity, respectively, in 
the notional elevation plane. Both planes are fixed relative to the spacecraft body 
frame. The parameters in Figure 17 are, again, assigned unit values. 

Common practice for implementing a slew maneuver to satisfy (3.10) is to 
use a so-called bang-bang or bang-off-bang trajectory where the maximum 
allowable force is applied over a short period of time to accelerate each movable 
joint, then the same magnitude of force is applied in the opposite direction to 
decelerate the links to their desired positions. This approach leads to the idea of 
computed torque control, in which required forces are generated from knowledge 
of the current state and its comparison to a desired profile. The control system 
uses feedback about the current state to determine the amount of control to 
apply [4], By this method, the spacecraft can maneuver along the prescribed 
path. A control law for such a system takes the form given in (3.11), which is a 
simple PD controller. 

f = K P { x d- x ) + K d{xd-x) ■ ( 3- 11 ) 

In (3.11), K p and K d are gains that can be programmed to achieve desired 

characteristics of the system. This method most likely produces a suboptimal 
solution but is easy to implement. Figure 18 shows the results of the maneuver 
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implementation using computed torque control for a rigid-body model; it also 
shows the behavior of the system after completion of the slew (indicated by the 
circles). Gains for the system are based on a one-second settling time with 
critical damping « = 1) to minimize overshoot. The control law in (3.11) is applied 
to both axes (azimuth and elevation) independently to achieve the desired final 
conditions in (3.10) within the desired time of 38.1 seconds. 
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Figure 18. Rigid-body slew maneuver using computed torque control: (a) 

position; (b) rate; (c) acceleration. 

To model the flexible-body motion of the system, the responses of the 
transfer function matrix in (3.2) can be simulated by applying the forces produced 
from the control law in (3.II). 2 Recall that X(s) = H(s)F(s), or 

X.(s) = H.,(5)F,(s) for j = 1,2,3 and k = 1,2,3 . Therefore, 


2 MATLAB's “Isim” command is very useful for this. 
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(3.12) 


JCl (5) = H 11 (5)F I (5) + H 12 (5)F 2 (s) + H 1 ,(5)F3(s) 

X 2 (s) = H2i(s)F|(s) + H22(5)F 2 (5) + H23(s) F 3(s) 
X 3 (s)='R 31 (s)F 1 ( S )+H 32 ( S )F 2 (s)+H 33 (s)F 3 ( S ) 

We are only concerned at this point with x 2 and x 3 (the control law will be set to 
maintain x 1 =0). Figure 19 shows the results. 



Figure 19. Flexible motion due to slew in Figure 18. 


Note how the peak amplitudes of the vibrational components coincide with the 
changes in acceleration in Figure 18, though the vibrations occur on a much 
smaller scale than the gross motion of the system. In addition, there is some 
flexible spillover. In other words, even though the gross motion is fixed beyond 
t f = 38.1, the system continues to oscillate in accordance with (3.12). 

Suppose now that it is desired to constrain vibrations to within 0.1 position 
units for both links. A difficulty in using computed torque control arises because 
one must set the appropriate limits on force and velocity to ensure this constraint 
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can be met. Finding these limits requires “tweaking” the specifications of the 
system until the desired effect is achieved. However, the total slew time can 
increase, which may not be desirable. For example, an increase in slew time of 
three seconds to 41.1 seconds reduces vibrations to below the constraint (see 
Figure 20). Figure 21 compares the acceleration commands for each case; the 
controller applies less torque during the slower maneuver, and this is the 
mechanism used to reduce the magnitudes of the oscillations. Hence, using a 
conventional computed torque strategy, the only way to control flexible effects is 
by manipulating the acceleration of the joints. This is why a flexible antenna must 
be slewed slowly to avoid such disturbances. 



Figure 20. Reducing the flexible effects by increasing the slew time. 


60 










Figure 21. Comparison of acceleration for different slew maneuver times. 

E. SUMMARY 

In this chapter, the mass-spring-damper model was expanded to include a 
third degree of freedom. This system is analogous to a spacecraft body with an 
antenna mounted on a double-axis gimbal. Conventional control techniques were 
then used to perform a notional antenna slew to illustrate the resulting flexible 
effects. Finally, it was shown that when using conventional control flexible effects 
can only be reduced by increasing slew time. There is no analytical process for 
determining the time increase that would properly constrain vibrations, so the 
process can be tedious. This simple exercise illustrates the state of practice; 
meeting constraints on vibrations is handled by operating the system more 
conservatively. However, this limits system throughput since slew time is wasted 
time. In the next chapter, optimal control is used to develop slew maneuvers that 
alleviate these issues. 
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IV. OPTIMAL CONTROL OF A 3-DOF SYSTEM 


It is desired to minimize slew time in order to maximize throughput. In 
addition, vibrational effects due to force inputs should be limited to their existing 
specifications. That is, moving faster should not perturb the system any more 
than the preset solution. Rather than use conventional control to minimize flexible 
effects, and potentially increase slew time, it is possible to determine an alternate 
control to minimize both flexibility and slew time. This can be done by 
manipulating the system through application of force in such a way that the 
flexible properties of the system may be used to an advantage. Finding such a 
solution can be done by formulating and solving an optimal control problem. 

Following [6], the problem must be redefined to satisfy the necessary 
conditions to achieve optimal control of the system. First, the dynamics of the 
system must be determined in state-space form. The correct definition of the 
state and control variables is crucial to a successful optimal control solution. 
Once the states and controls are properly defined, we can proceed to redefine 
our problem statement. 

To construct an optimal control problem, a cost functional /(x(*),u(»)) to 

be minimized should be identified. The cost, a function of the state vector x and 
control vector u in the system, consists of an endpoint cost (expressed in terms 
of the final state values) and/or a running cost (calculated for each time step). 
Thus, the cost functional can be expressed as 

^(xW.u(')) = ^(x(^)) + JV(x(0,u(0^ ■ (4.1) 

For minimum-time solutions, the cost functional is simply J = t f , but variations of 
this simple objective function are described below. 
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The next step is to define any boundary conditions for the states, controls, 
and time. Any constraints on states or controls, (e.g., |/ max | = 1.5), can also be 
applied. The redefined problem statement could be written as 


Minimize 
Subject to 


x = f(x(f),u(f)) 



*i(°) * 2 (°) 



t 0 =0 

e ( x M) = 0 

x L <x i <x u , i = 1 , 2 ,... n 
u L < u ( . < u y , i = 1 , 2 ,... m 


(4.2) 


To solve (4.2), the Hamiltonian must be derived. The Hamiltonian is a function of 
the states, costates, and controls, where the vector of costates is related to the 
dynamics of the system [6], The Hamiltonian is defined as 

H(x,X,u) = F(x,u) + A r f(x,u) (4.3) 


where X are the costates. 


The necessary conditions to achieve an optimal control solution can now 
be determined [6], The first necessary condition is the Hamiltonian minimization 
condition, written as 


r)H 

3u 


= 0 . 


(4.4) 


The second necessary condition is the satisfaction of the adjoint equation(s) 


i(t) = 


dH 

dx 


(4.5) 


The transversality condition relates the final values of the costates to the 
Endpoint Lagrangian and comprises the third necessary condition, written as 



(4.6) 
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where the Endpoint Lagrangian £’(v,x(r / )) = £’(x(r / )) + v r e(x(r / )) [6], 

The MATLAB software package DIDO automatically determines an 
optimal control solution that satisfies the above conditions [6], DIDO is used to 
solve the optimal control problems posed through the rest of this thesis. 

A. STATE-SPACE FORM OF COUPLED SYSTEM 

To construct the optimal control problem, the system in (3.9) must be 
transformed to state-space form. Assuming unit values for mass, stiffness, and 
damping, let y e M 3 * 1 be the total displacement of each mass such that 


Ti 


Xj + X j 

y 2 

= 

X2 H - X2 

y 3 


X3 + ^3 


Let z e M 6j1 be a state vector such that 


z = 


y 

y 


Xj + X j 

%2 X 2 

X3 + X 3 
Xj + X^ 

%2 X-2 

X 3 + X-2 


(4.7) 


(4.8) 


Let u e M 3j1 be the control vector u = 


fi fi f 3 


I T 


(i.e., the vector of forces 


acting on each mass). The states represent the total motion described by the 
sum of the rigid-body component and the flexible-body component, as described 
by (3.9). The generic state-space representation for a rigid-body slew is 


X 

X 


0 1 
0 0 




(4.9) 
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where X is the displacement, X = v is the velocity, X = a is the acceleration, m 
is the mass, and f is the control force. This state-space form, adapted for the 
3-DOF system, is 


*1 

Z 2 

*3 

X 2 

*3 


0 0 0 1 0 0 
0 0 0 0 1 0 
0 0 0 0 0 1 
0 0 0 0 0 0 
0 0 0 0 0 0 
0 0 0 0 0 0 


0 0 




(4.10) 


where u = ^ 0 f 2 / 3 J since the restoring force is due to the control law on 

mass 1 comprising k x and c, (which will be accounted for as part of the flexible 
dynamics). The flexible-body component in state-space form becomes 





(4.11) 
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Developing the state-space representation using (4.8) as the state vector is 
grueling. One cannot simply add the rigid- and flexible-body state-space 
equations together to produce a form that is conducive to modal analysis. The 
difficulty lies in that the flexible motion equations contain no X terms or their time 
derivatives. For example, the equation for z 4 is 

z 4 = x x + X, = 2/j - 2x 1 + x 2 - 2i, + x 2 . (4.12) 

Instead, separate the rigid-body and flexible-body components and define a new 
state vector. Let y e R 12vl be the state vector 

~\ T 

y= X, X 2 X 3 X, X, X 3 x l x 2 x 3 i, x 2 x 3 (4.13) 

and 

-i T 

y= X, X 2 X 3 X, X 2 X 3 ij x 2 x 3 x x x 2 x 3 . (4.14) 

The state-space representation, with unit values substituted, is therefore 

0001000 0 0 0 0 0 1 [o 0 0 

000010000000 000 

000001000000 000 

000000000000 100 

000000000000 010 

0000000 0 0 0 0 0 001 (4 15) 

0000000 0 0 1 0 0 000 

0000000000 1 0 000 

00000000000 1 000 

000000 -2 1 0 -2 1 0 100 

0000001 -2 1 1 -2 1 010 

0000000 1 -1 0 1 -lJ [ o 0 1 

Equation (4.15) provides the dynamics required for the optimal control problem 
definition. Recall that (4.15) represents coupled motion in physical space. The 
time-optimal slew given the conditions in (3.10) and dynamics in (4.15) can be 
determined by solving the following problem: 
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Minimize / = t f 

Subject to v, = y 4 

y 2 = y 5 
y 3 = y 6 

y 4 = m, 

y$ = u i 

y 6 = u 3 

y 2 = Jio 
^8 = 1 

j 9 = _Vi2 . (4.16) 

-Vio = -2^7 + y 8 - 2 -VlO + Vu + “l 

yi ,=y 7 -2y 8 + y 9 + y 10 - 2y,, + y 12 + u 2 
yn=y*-y9 + yu-yn+ u 3 
(y 2 {h),y i {t 0 ),y 5 {t 0 ),y 6 {t 0 )) = (-lO,5,0,Q.l) 

(-V 2 ( ? /)^3(^/)^5 (^/)) = ( 15 -15,0.1,0) 

h) = 0 

|y,!< 0-75,/ = 4,5,6 
\u\ < 1.5,i = 1,2,3 

The dynamics in (4.16) are such that constraining the total motion of the 
system is simple. However, they do not readily provide modal characteristics. 
The dynamics are expressed in physical space, so the modes cannot be 
constrained directly. The dynamics must be uncoupled by transforming them to 
modal space, thereby enabling effective mode constraints. 

B. UNCOUPLING THE 3-DOF SYSTEM 

The system as constructed exhibits proportional damping. Therefore, the 
undamped mode shape vectors can be used to diagonalize the mass, damping, 
and stiffness matrices to produce the modal matrices, and thus the system can 
easily be uncoupled into three single-DOF systems of the form 

ih,.q r + c r q r + k r q r = .// , (4.17) 
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where r is the mode number. Recall from Chapter II that f'= V r f. The transfer 
functions in the Laplace domain can then be found easily. Assuming 

q(0) = q(0) = 0, 

r k F' 

Q r s 2 + ^rQ r s + ^rQ r =^r 
m r m r m r 

] _ (4.18) 

q ,.( 5 )= n r (s)¥' r (s) = — m ‘ _ f ; c . v ) 

s 2 + C ^s H-— 

m r m r 

These transfer functions represent displacement in modal coordinates due to the 
transformed forces, and they are much easier to derive and manipulate than the 
transfer function matrix in (3.2). It is important to distinguish the difference 
between an uncoupled transfer function in modal (principal) coordinates versus a 
coupled transfer function in physical coordinates. 

The uncoupled equations of motion from (4.17) are 

f' 

q x +0.1981<7, +0.1981^ = ——— 

f' 

q 2 + l.5550q 2 + l.5550q 2 = Jl . (4.19) 

1.8412 

f' 

q 3 +1.5550<7 3 +1.5550g 3 =--— 

2.8629 

These equations represent uncoupled modal masses as in Figure 22. The 
optimal control problem definition can now be redefined in terms of modal 
coordinates. Let the state vector be defined as 


y = 


X\ X 2 X 3 Xj X 2 X 3 q t q 2 q 3 q A q 5 q 6 


where 


^4 c h 



(h q 2 Qi 


T 

, and 


-]T 


x, X 2 X 3 X t X 2 X 3 q ! q 2 q 3 q 4 q 5 q 6 


(4.20) 


(4.21) 
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; m 
q 3 = ° 


Figure 22. Uncoupled 3-DOF mass-spring-damper system. 

-| T 

The control vector remains u= 0 f 2 f 3 . However, the dynamics equations 

-i T 

for the flexible component defined by % ••• y n must account for the 

transformation from physical coordinates to modal coordinates. For example, 
from (4.19) we can see that 

^ 1 = 9 ^ 9 _0 - 19814l_0 - 1981?1 ' (4 ' 22) 

Because f'= v u f + V 2 J 2 + V 3 J 3 = m ] (v u a t + V 2 i a 2 + V 3 l a 3 ) and q l = q 4 , the general 
form is 

q4 = ^ L (V u a 1 + V 21 a 2 + V n a 3 )--^q 4 -^q 1 . (4.23) 

m, m l m, 

Thus, the general form of the optimal control problem in modal coordinates can 
be defined as in (4.24). Note that there are no constraints or terminal boundary 
conditions placed on the flexible motion in the problem statement. The effects of 
constraining the flexibilities are discussed in a later section. 
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Minimize 
Subject to 


J = t f 

y, = y 4 
y 2 = y 5 
y 3 = y 6 

u. 

y 4 = — 

m i 

«, 

m 2 


m 3 

y i = y w 

% = -V, 1 


.(4.24) 


y 9 = y 12 

1 C k 

y W = — ( y .l«l + y 21«2 + V 31« 3 ) - “Jv 

m, m j ra, 

1 c k 

j'l. = — + ^22«2 + y 32« 3 )- y 8 

m 2 m 2 m 2 

1 c k 

?12 = — ( y !3 M l + y 23«2 + ^33 M 3 )- TT->'l2 - 

m 3 ra 3 m 3 

(} ; i( ? o)^ 2 ( ? o)>> , 3( ? o)>3 ; 4( ? o)»>’ 5 ( ? o)>>'6( ? o)) = ( 0 -10,5,0,0,0.1) 
(yi(^)>y 2 (^).y3(^).y4(^).y s (^).y 6 (^))=(o,i5,-i5,o,o.i,o) 

(3 , 7( f o)’3 , 8 ( f o)’3 , 9( f o)’3'lo( r o)’3 , ll( f o)’3 ; i2( f o)) = ( 0 - 0 - 0 - 0 - 0 ’ 0 ) 

|y ; | < 0.75, / = 4,5,6 
|u,.| <1.5, i = 1,2,3 


C. SLEW OPTIMIZATION WITH NO CONSTRAINTS ON VIBRATIONS 

The MATLAB package DIDO [6] was used for optimizing and validating 
the minimum-time maneuver. The minimum-time maneuver is presented in 
Figure 23. For both axes, the maximum acceleration (instead of some fraction of 
it) is supplied initially and at the end of the maneuver for a much shorter time 
than in the computed torque case. It is worth noting that y 3 does not follow the 
exact same trajectory as in the computed torque slew in Figure 18. The 
acceleration is not supplied in a “bang-off-bang” method but rather oscillates 
slightly. Because y 2 (r / )-y 2 (0)>_y 3 (r / )-_y 3 (0), the azimuth axis becomes the 


limiting axis. As a result, the elevation axis has additional freedom in the path it 
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takes to the desired endpoint. This is why the acceleration applied to that axis is 
not uniform. As an analogy, picture two baseballs being thrown simultaneously 
from the same initial position with the same constraint on maximum velocity, but 
with the first baseball’s target 10 feet further than the second’s. The shortest 
travel time for the first baseball dictates the shortest travel time total for both 
baseballs together. The second baseball’s velocity does not need to match the 
first baseball’s velocity to reach its target at the same time. For the optimal 
control problem in (4.24), that time is 33.8 seconds, a savings of 4.3 seconds 
(12.7%) over the computed torque maneuver. 



Figure 23. Rigid-body results of optimal control problem defined in (4.24) 
(t f = 33.8): (a) position; (b) rate; (c) acceleration. 


Figure 24 shows the resulting flexible motion in both modal coordinates (in 
terms of q) and physical coordinates (in terms of x). Recall that x = Vq and that 
q r ^ x r . Rather, q l represents the motion of all masses with respect to mode 1, 
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while Xj is the physical displacement of mass 1 only, and so on. Figure 24 also 
illustrates the relationship x = Vq. For example, 

q(5) = [ -0.1165 -0.0231 0.0051 J , and the second row of the mode shape 

matrix is [ 1.8019 0.4450 -1.2470 ]. Thus, x 2 (5) = -0.2266, as the plot shows. 
It is also evident that mode 1 dominates the total motion of the system. 



Figure 24. Flexible-body results of optimal control problem (4.24): 
(a) modal coordinates; (b) physical coordinates. 


The added freedom in the path of y 3 (discussed previously) allows control 
designers to apply cost reduction measures to the system. For example, it is 
possible to find a minimum-time maneuver that also minimizes the energy 
expended to control the non-limiting axis. In this case, the problem definition in 
(4.24) can be altered to minimize the cost function 
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(4.25) 


f / 

J = 200 J u 2 dt + t f 

0 

subject to the same dynamics and boundary conditions/constraints. The integral 
term in the cost function represents a running cost, in that it is calculated and 
summed at each time step, and it is weighted 200 times against the endpoint 
cost of final time. By doing this, the control (energy) is reduced while also 
minimizing time. The new angular rates and accelerations for the revised optimal 
slew are shown in Figure 25. The displacement of y 3 now more closely matches 
the computed torque slew while completing the maneuver in the same time as 
Figure 23. Due to (4.25), the maneuver realizes additional savings in energy 
expended. 



0 5 10 15 20 25 30 35 

(C) 

Time (sec) 


Figure 25. Rigid-body results of minimum-time slew with running cost to 
minimize control effort [t f = 33.8) : (a) position; (b) rate; 

(c) acceleration. 
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The new acceleration profile ostensibly affects vibrations, shown in Figure 26 and 
Figure 27. In Figure 26, mode 1 again clearly dominates the total motion of the 
system but in a more pronounced way as the spacecraft minimizes its control 
input. Modes 2 and 3 show less oscillation than before (compare Figure 24 to 
Figure 26). 



Figure 26. Flexible-body results of optimal control problem with running 

cost in modal coordinates. 

Figure 27 shows a comparison of each mass’s displacement (magnitudes only) 
for each optimal control solution (minimum-time and minimum-time plus minimum 
control effort). The addition of a running cost generally reduces the magnitude of 
displacement across all three masses, which follows as the applied control is 
minimized. Applying a running cost to minimize applied control can be a feasible 
method of reducing vibrations through a maneuver and is one of the tools 
available to designers using the optimal control approach. Note that this option is 
not possible using conventional control concepts. 
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The plots of the modal coordinate q illustrate an interesting phenomenon. 
Recall from (3.6) that the damped frequencies for each mode are 

" 0.4339 0.9749 0.7818 ] rad/s, or [ 0.0691 0.1552 0.1244 ] Hz. The 

natural tendency for any system is to oscillate at its natural frequency [7], The 
same is true for this system. In Figure 26, for example, the period of q 2 
decreases over time before the last acceleration pulse is applied, eventually 
reaching about 6.5 seconds, or roughly 0.154 Hz. 


Xi X 2 



.Minimum time only 

-With running cost 



Figure 27. Comparison of flexible-body results in physical space for 
endpoint cost vs. endpoint plus running costs (magnitudes only). 


The implementation of the spacecraft’s control system also becomes 
apparent. The displacement of stays at or very near to zero (on a scale of 

10 10 , see Figure 23) throughout the rigid-body simulation but its flexible-body 
counterpart does not. This vibrational motion is analogous to the control system 
counteracting the effects of the motion of the two axes to keep the spacecraft 
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stable in its position. As we constrain the control input, and thus reduce the 
applied force on y 2 and y 3 , it follows that the relative magnitude of q l decreases 
since mode 1 is closely mapped to the rigid-body motion of the system (as 
previously discussed). By reducing the force used to slew masses 2 and 3, the 
spacecraft’s control system does not need to work as hard to counteract the 
displacements of the other axes. 

Due to the excitation of the flexible modes, the spacecraft and its 
appendages will continue to vibrate after the slew has completed (see Figure 25). 
Figure 25 shows the behavior of the vibrations after the final conditions are met, 
when t > t f and the payload (antenna) must be operated. In Figure 28, the 

physical properties of each mode take over as the slew control system 
deactivates. Modes 2 and 3 dampen quickly while mode 1 takes longer to settle. 



Figure 28. Flexible spillover for time-optimal slew maneuver: (a) modal 

coordinates; (b) physical coordinates. 


Further, referring to Figure 28, as mode 1 settles, its frequency approaches its 
damped natural frequency. The bodies in each axis continue to physically 
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oscillate about equilibrium but at the same frequency (dominated by mode 1). 
From this result, it is apparent that performing the slew maneuver excites the 
base body in a way that perturbs the motion of the entire system. If the 
perturbation is large enough, it may be the case that the operation following the 
slew, (e.g., communicating over an antenna link), could be disrupted. Therefore, 
it is desirable to try to limit this flexible motion. This will be the subject of the 
sections that follow. 

D. SETTING TERMINAL CONDITIONS ON THE FLEXIBLE MOTION 

Suppose it is desired to extinguish all vibrations at the end of the 
maneuver. The appropriate terminal flex constraints added to the problem 
statement in (4.24) are 

(<5r ] (r / ), ?2 (r / ), ?3 (f / ),^4(f / ),^(f / ), ?6 (r/)) = (0,0,0,0,0,0) - (4.26) 

One would expect that more control or more time is required to cease flexible 
motion and meet the desired end conditions for the maneuver. The minimum¬ 
time solution is shown in Figure 29 and confirms this assertion. The maneuver 
now takes 35.14 seconds, an increase of 1.34 seconds (4%), as more time is 
required to null the vibrations. Similar behavior occurs in the elevation axis as 
discussed previously since the spacecraft must travel a greater distance in the 
azimuthal direction. Note that in Figure 29 that the trajectories are plotted at the 
solution nodes. They would be smoothed by the correct interpolation prior to 
implementation on a real system. 


78 




Figure 29. Optimized maneuver with zero-flex conditions at endstate 
[t f = 35.14): (a) position; (b) rate; (c) acceleration. 


Figure 30 plots the maneuver in modal and physical coordinates. Note the 
impulse on all three axes at ~23 seconds; this coincides with the acceleration 
impulse along the elevation axis in Figure 29(c) and provides an illustrative 
example of how the three axes are coupled. Force is applied on only one mass, 
yet this induces motion on the other masses. In Figure 30(a), mode 1 again 
dominates the total motion. 
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Figure 30. Flexible spillover with terminal flex constraint: (a) modal 

space; (b) physical space. 


A running cost is now added to minimize control. Applying the cost 
function in (4.25), the optimized maneuver (Figure 31) takes 35.26 seconds, only 
marginally slower than in the previous case (0.12 seconds, 0.36%). This shows 
that it is possible to reshape the trajectory of the non-dominant slew axis in order 
to reduce the control effort. Note that the acceleration profile for the azimuth axis 
(y 2 ) is essentially the same as before. However, the reduction of the control 
effort on the elevation axis has an advantageous impact on the flexible effects. 
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Time (sec) 


Figure 31. Maneuver optimized for time and control effort [t f = 35.26): 

(a) position; (b) rate; (c) acceleration. 


Figure 32 illustrates the flexible components of the minimum time plus 
minimum effort slew. The total magnitude of flexible-body displacement has been 
significantly reduced. Variable x 3 alone experiences roughly a 50% decrease in 
peak displacement. With the negligible increase in maneuver time and decrease 
in control and vibrations, this maneuver proves to be more desirable from the 
point of view of system operation since the load on the base body controller is 
reduced. 
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Figure 32. Flexible effects from the maneuver in Figure 31: (a) modal 

space; (b) physical space. 


An interesting phenomenon occurs in these two maneuvers where the 
system uses the additional freedom in the motion of the elevation axis to meet 
the desired endstate while maintaining steady acceleration on the azimuth axis. 
The controller excites the modes of the system, by applying force only on m 3 , in 
such a way that its modal properties work to reduce the effects of flexibility in the 
system. This non-intuitive aspect is only possible by using optimal control to 
design the maneuver trajectories. 

To further illustrate the point, assume that the initial conditions are 
changed to 

(Ti (^o)’T 2 ( ? o)’T 3 ( ? o)’T 4 ( ? o)’T 5 ( ? o)’T 6 ( ? o)) = (0,-to,to, 0 , 0 , 0 . 1 ) (4.27) 

so that both axes now must move the same distance, such that 
.^2 ( r /) — (°) = - v 3 ( r /) — -Ts(°) ’ thus removing freedom previously present in the 

elevation axis. The results are telling; the minimum-time maneuver (no running 
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cost) is 0.89 seconds (2.5%) slower (see Figure 33), and the controller must now 
torque both axes before the end of the maneuver in order to suppress the 
flexibilities. A similar effect is observed for the minimum time and effort problem 
(see Figure 35). The vibration profiles for each maneuver (see Figure 34 and 
Figure 36) are similar, further emphasizing the role a non-dominant axis can play 
in accommodating flexible effects. 



(b) 



Figure 33. Minimum-time maneuver for 

(t f = 36.03): (a) position; (b) rate; (c) acceleration. 
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Flexible displacement 



(b) 

Time (sec) 


Figure 34. Flexible displacement from minimum-time maneuver of 
Figure 33: (a) modal space; (b) physical space. 
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(b) 



Figure 35. Maneuver optimized for time and control effort with 

y 2 (^/)-y2(0) = ^ 3 (^/)-y 3 (0) = 36.13): (a) position; (b) rate; 

(c) acceleration. 
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Figure 36. Flexible displacement from maneuver of Figure 35: (a) modal 

space; (b) physical space. 


E. BOUNDS ON FLEXIBILITY DURING SLEW 

A less stringent case than designing for zero spillover at the end of a 
maneuver is to set a bound on the flexible displacement during the maneuver. In 
other words, maintain flexible motion within a desired range. Perhaps a designer 
is willing to trade flexible effects for slew time or control effort. One cannot simply 
set nonzero bounds on q(r) and expect them to translate into the physical world. 

As discussed previously, the transformation x = Vq must be applied. For 
example, suppose it is desired to bound vibrations in physical space for all axes 
to the interval [-0.3, 0.3] degrees with no constraint on x. The transformation for 
the upper and lower bounds yield, for the example problem of this chapter, 

q uppe r=[ 0.1629 0.1048 0.0323 J (4.28) 
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(4.29) 


q w = [ -0.1629 -0.1048 -0.0323 J . 

However, the simple box bounds given in (4.28) and (4.29) do not properly 
constrain q, resulting in potential violations of the constraint on x (see 
Figure 37). 



Time (sec) 

Figure 37. Flexible displacement for a slew with constraints on x applied 
via x = Vq: (a) modal constraints satisfied; (b) physical 
constraints violated. 

The constraint violation in Figure 37 occurs because each modal coordinate 
depends on the values of all three physical coordinates (and vice-versa when 
transforming the other direction). For a simple explanation, consider a 2-DOF 
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system with eigenvector matrix V = 


. A plot of the box bounds on x is 


a b 
c d 

shown in Figure 38 along with corresponding bounds on q = V 'x. 


“Naive” bounds 



Figure 38. Example constraints on x transformed to q via q = V 'x. 


The area in the figure outlined in solid black represents all of the values of q } and 
q 2 in modal space that satisfy the constraints in physical space. However, the 
transformation of the bounds on x directly into modal space produces the large 
rectangle in dashed red (naive bounds). In this case, there are points in modal 
space that satisfy the constraints on x but not on q. Plots of q transformed to x 
would show the same behavior. As seen in Figure 37, a transformation of bounds 
in physical space to modal space and back results in violation of the constraints. 

One approach to overcome this is illustrated by the small rectangle in 
dashed red in Figure 38. Instead of constraining x 1 and x 2 simultaneously, 
perform a transformation on each individually and use the minimum absolute 


88 




















values of q to form a conservative constraint boundary as in Figure 38. For the 
example 3-DOF system, the modal constraints {q x } for each axis are 

q Xi = V“' [ 0.3 0 0 J=[ 0.0323 0.1629 0.1048 ]' 

q X2 = [ 0 0.3 0 J=[ 0.0582 0.0725 -0.1307 ] ? ■ (4.30) 

q Xi = V" 1 [ 0 0 0.3 J=[ 0.0725 -0.1307 0.0582 J 

The resulting constraints on each modal coordinate are 

q„,-,«=±[ 0.0323 0.0725 0.0582 ]' (4.31) 

This approach creates conservative bounds in modal space but ensures that the 
maneuver stays within desired bounds in the physical space. Interestingly, when 
optimizing the maneuver for both time and control effort, the maneuver looks 
similar to the ones generated using bang-off-bang control (see Figure 39). The 
maneuver time is 34.50 seconds. The induced vibrations (see Figure 40) are 
observed to ride the constraints indicating the conservatism inherent to the 
maneuver. 


89 




0 5 10 15 20 25 30 35 

(C) 

Time (sec) 

Figure 39. Maneuver optimized for time and control effort with 
conservative constraints on q [t f = 34.5o): (a) position; (b) rate; 

(c) acceleration. 
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Figure 40. Flexible displacement for the maneuver of Figure 39: 
(a) modal space; (b) physical space. 


As an alternative, a path constraint can be used to implement equations 
that describe the true modal constraints shown in Figure 38. In DIDO, the 
constraints are applied in physical space by transformation of the current modal 
state, as 


p = V 


yi y% y 9 



(4.32) 


Using (4.32), the constraints in physical space are applied to the transformation 
of the modal coordinates at each node. The optimal maneuver with the path 
constraints is shown in Figure 41 and Figure 42. Note that the execution time is 
close to that of the previous case, and additional torque must be applied to each 
axis at certain instants in time to avoid violating the constraint placed on the 
amplitude of the vibrations (compare Figure 39(c) with Figure 41(c)). Recall that 
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the bound in physical space vibrations is -0.3 <x <0.3; the path constraint 
properly enforced this bound as shown in Figure 42(b). 



(b) 



Figure 41. Maneuver optimized for time and control effort, using a path 
constraint to suppress vibrations (t f = 34.54): (a) position; 

(b) rate; (c) acceleration. 
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Figure 42. Flexible displacement for maneuver in Figure 41: (a) modal 

space; (b) physical space. 


These slew maneuvers are compared to the ones from Chapter III in 
Table 1. Not only is it possible to constrain flexible motion, but slew times 
decrease when using optimal control rather than increase when using computed 
torque control. The table highlights the fact that optimal control provides the 
designer with multiple options, each reducing slew times compared to traditional 
control methods. 
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Table 1. Comparison of slew maneuvers using computed torque 

control and optimal control. 



Computed torque 
control 

Optimal control 

Time 

savings 

h 

(sec) 

Max amplitude 
(deg) 

h 

(sec) 

Max amplitude 
(deg) 

Relative to 
computed 
torque 

During 

slew 

Spillover 

During 

slew 

Spillover 

Unconstrained flexibility 

38.1 

0.143 

0.118 

33.8 

0.392 

0.266 

11.3% 

Zero terminal flex 

Cannot execute 

35.14 

1.04 

0 

7.7% 

Zero 

terminal 

flex 

Conservative 

bounds 

Cannot execute 

34.50 

0.15 

0 

9.5% 

Path 

constraint 

Cannot execute 

34.54 

0.30 

0 

9.3% 


F. SUMMARY 

In this chapter, the antenna slew problem was reformulated in order to 
apply optimal control theory to determine the control effort required to minimize 
slew time while simultaneously managing structural flexibility. It was shown that it 
is possible to apply torque in such a manner that the flexible properties of the 
system could be used advantageously to control vibrational effects. It was also 
shown that optimal control could be utilized to minimize the control effort without 
significantly increasing slew time or flexible motion. This effect is particularly 
useful in instances where one link has a greater distance to travel than the other. 

This chapter additionally illustrated the value of working in modal space 
versus physical space. Modal space readily produces the natural properties of 
the system by mode and allows control system designers to understand the 
responses of a flexible system to force inputs. The equations of motion in modal 
space represent single-DOF systems that uncouple the physical space 
equations, making analysis and control of the system dynamics much easier. 

The challenges of defining the constraints on vibrations in physical and 
modal spaces were discussed. The key difficulty is preserving the relationship 
x = Vq. Conservative bounds in modal space ensure that physical constraints 
are met, but the system must exert more effort than is required to stay within 
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those bounds. Naive bounds on system behavior in modal space may cause 
violation of the bounds in physical space. Therefore, a path constraint that uses 
the relationship x = Vq must be used to ensure the correct physical constraints 
are met. 

The optimal control solutions presented in this chapter greatly reduced 
flexible effects while simultaneously decreasing slew time. While computed 
torque control from Chapter III required constant tweaking of the system to 
achieve desired results, the optimal control method allows for easy constraint 
definition as long as the problem is posed correctly. The optimal control approach 
will be applied to a more realistic antenna system in the next chapter. 
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V. OPTIMAL CONTROL OF A DOUBLE- 
GIMBAL ANTENNA WITH FLEXIBLE JOINTS 


This chapter introduces a double-gimbal mechanism (DGM) to add to the 
system from Chapter IV to provide a more realistic model of an actual gimbaled 
spacecraft antenna. The equations of motion of this system are significantly more 
complex than those of the simple models in previous chapters. The rotational 
dynamics introduce additional nonlinear terms that are dependent on the current 
state of the system. Building on the developments of the previous chapter, this 
chapter compares the results from computed torque control and optimal control 
of the nonlinear system. The goal of this chapter is to show that modal analysis 
and optimal control can be applied effectively to nonlinear systems. 

A. UPDATING THE MODEL 

Figure 43 shows a DGM that is now applied to the model. The antenna is 
attached to the spacecraft body through this two-axis gimbal. 



Figure 43. Diagram of two-axis gimbal (body 1 represents azimuth, body 

2 represents elevation), from [5], 
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The equations of motion for this rigid-body system are derived in [5] as 


/ i^ + / 2x r sin 2 0 2 +/ 2:: cos 2 0 2 

0 

0 . 



0 

1 

Si 

CN 

1 

CN 

_ i 


1 

to 

1 _ 


2 ( / 2«- / 2 ZZ )^Asm0 2 cos0 2 

c°s0 2 


(5.1) 


In (5.1), Oy and 0 2 represent the angular displacement of the azimuth and 
elevation, respectively. This model does not account for the motion of the 
spacecraft body and ignores flexibility. For example, the motion of the azimuth 
joint coincides exactly with the appendage connected to that joint. 

The model of (5.1) is now extended to include flexible effects. Figure 44 
shows a schematic of a flexible joint. The spring with stiffness constant k t 
represents joint flexibility and allows elastic displacement between the rotor 
(mounted on the spacecraft) and the rotary link. The dashpot with damping 
constant c ; represents the torsional damping due to rotational friction, structural 
properties, and external mechanisms (in other words, anything that dissipates 
energy during motion instead of storing it). 



Figure 44. Representation of an elastic joint, after [14]. 


98 










In Figure 44, y 2i is the angular displacement of axis i, and y 2M is the angular 
displacement of the actuator. These displacements are related to the link and 
rotor angles as 


y 2i =0 i ,i = l,...,n 


Yu-i =-—■n 


where n is the number of axes in the system (in this case, n = 2) and G ( is the 
gear ratio of the actuator. For the purposes of this model, it is assumed that the 
system consists of a direct drive between the actuator and shaft, so there are no 
gears between the actuator and rotor. As a result, from (5.2), y 2M =0 ; -. The 

elastic displacement in the joint, therefore, is y 2i -72 M- 

r ~\ T r ~\ T 

Let k t = k, c t = c, y,=[ y 2 y 4 J =[ d l 9 2 J , and 

-| T r -| T 

y 2 = 7, y 3 = 0i 02 ■ The inertia tensor of the rotor is given as 

(l, OM ,),=diag[ l, xx l iyy l izz ], and D(y,) gives the configuration-dependent 

inertia of the rigid-body. Using these definitions, the kinetic energy for the system 
may be computed as 

KE = ^y\D{y l )y, + ^y T 2 i roto j 2 (5.3) 

where moment of inertia J rolor = diag Gf(l rotor ) la ••• G;(l rotor ) nzz [14]. The 

potential energy, a function of the angular displacements of the axes and the 
elastic displacement, is 

PE = P 1 (y i ) + P 2 (y 1 - y 2 ) (5.4) 

where function P l is the potential energy of the rigid-body system. Function P 2 
depends on the elastic displacement and is defined by 
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(5.5) 


1 T 

p 2 = 2 k {y t ~y2) (yi-y 2 ) ■ 

The Lagrangian may be constructed as 

L ~ KE — PE = |(yfZ)(y,)yi + y 2 J rotor y 2 - 2 P x (y t ) - k( y, - y 2 f (y, - y 2 )) (5.6) 

Damping results in energy dissipation rather than energy storage. A typical 
model of this is Rayleigh dissipation given as [15] 

^^■(y.-y,) 2 (5.7) 

In the Lagrangian of (5.6), dissipation may be added to the total energy, and thus 
the Euler-Lagrange equation becomes 


' OL] I dfffa - 

v 3 y,J By, 5 y 


The generic equations of motion may now be derived from (5.8) as 

^yi + f^y,-^-yf|^-yi + |^-l+c(y 1 -y 2 )+fc(y 1 -y 2 ) = 0 (5.9) 

l 2 dy, dy ,) 

J, otor y 2 -c(y 1 -y 2 )-^(y 1 -y 2 ) = u . (5.10) 

Equation (5.9) governs the motion of the rigid links, while (5.10) governs the 
motion of the rotors. The second term of (5.9) in parentheses represents the 
Coriolis, gravitational, and centripetal forces acting on the system. The new 
r i T 

control vector u= t, t 2 now acts on (5.10); its effects are felt in (5.9) via 
transmission through the flexible joint. 
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The new equations for the flexible DGM are thus 


A,, + 4v>n 2 0 2 +4;Cos 2 0 2 

0 

0. 

0 

hyy 
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<N 

_1 


2 ( / 2x,- / 2,,)^1^2 Sin0 2 COS0 2 

+ c 

01-01 

+ k 

1 

1 Jfc 
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(ha-hnWl sin0 2 cos0 2 


d 2 — 0 2 


^2 T2 


(5.11) 
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where the link torques in (5.1) are replaced by the torque transmitted via terms 

-c(0i-0,)-£(fli-0,) and -c(0 2 - 0 2 )-&(0 2 - 0 2 ). 


B. SIMULATION OF THE COUPLED MOTION 

Following from the example from [5], let / 1 = diag[ 4 4 6 ], 

/ 2 = diag [ 224 ], (l rotor ). = diag[ 111], and c = k = 1 with appropriate 
units. The equations of motion may thus be rewritten as 


6 + 2 sin 2 0 2 + 4 cos 2 0 2 0 

0, 

+ 

-40,0 2 sin0 2 cos0 2 

1 

(N 

O 

_l 

1 - 

to 

1_ 


20 2 sin 0, cos 0, 


(5.13) 
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A simple PD controller with closed loop feedback was used to simulate a 
conventional slew of the flexible system. The following control law was used. 

r = K l ,{e J -e)+Kja J -a) (5.15) 
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where O d and co d are the desired angular displacements and angular rates for 
each axis. Given a natural frequency co n = 0.6 rad/s and damping ratio £ = 1 to 
prevent overshoot, the gains may be computed in real-time as the inertia 
changes for the system (this is a version of computed torque control). 


Z P =< 


/ u z + 4«sin 2 0 2 +/ 2zz cos 2 0 2 


l 2yy 




r ^ + h^ s ^S 2 +I 2 a cos-B 2 


l 2yy 


(5.16) 


The slew was implemented with the following boundary conditions and 
constraints: 

(#i (*o )A (^o )»®2 ( ? o )) = (40°,20°,0,0) 
(0 1 (^).0 2 (^ / ).o 1 (^).© 2 (^)) = (O,4O o ,O,O) . (5.17) 

(r 0 ,r / ) = (0,20)s 

The damping and stiffness values affect how much (or little) actuation 
(control torque) is needed to slew each axis to the correct position. If the stiffness 
value k °°, then the elastic displacement 6 t -0 ( —>0. The joint now behaves as 
a rigid body. A slew for the rigid-body system is shown in Figure 45. Note the 
similarities in the torque profiles for the azimuth axis for this system and the rigid- 
body systems studied earlier in Chapter III. More effort is required to control the 
azimuth axis versus the elevation axis because of the time-dependent nonlinear 
inertia terms. 
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Figure 45. Computed torque control of flexible DGM with k » 0 : 

(a) angular displacement; (b) angular rates; (c) control torques. 

Figure 46 shows the angular displacements and rates of each axis and the 
torques applied to the rotors for k = 1, a much more flexible system. The inertia 
term for the azimuth axis changes as the elevation angle changes, affecting the 
gains and torques needed to control the system. Note how the controller 
responds to the changing inertia in the azimuth axis (see Figure 46(c)). As noted 
previously, it takes much more effort to control the azimuth axis than the 
elevation axis, but the torque profiles follow the same basic shapes as in 
Figure 45. Due to the elasticity in the system, the controller must continue to 
work even after t f in order to damp out the vibrational effects. The large flexible 

displacements after the slew (post-maneuver spillover) could present a problem 
for antenna operations. 
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Figure 46. Computed torque control of flexible DGM with k = 1: 

(a) angular displacement; (b) angular rates; (c) control torques. 


Figure 47 illustrates the relationship between the axis angle 9 and the 
rotor angle (j ). The elasticity in the joint causes the axis to lag the rotor as the 
actuator torques the rotor. The rotor must displace a certain amount before the 
axis follows suit. As stiffness increases, the time histories of each angle 
approach the same values as the lag decreases. 
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Figure 47. Time history of rotor, (j ), and link, 9, motion for a conventional 
slew with k = 1: (a) angular displacements; (b) angular rates. 

Note that in using (5.15) it is assumed that the actual position of the 
gimbal links can be measured. An alternate control law could be implemented if 
the positions must be inferred at the rotor. In this case, the torque law would be 

* = K p {O d -<l>) + K d (G>d-i>) ■ ( 5 - 18 ) 

Due to the elasticity and coupling in the system, the slew maneuver using (5.18) 
takes significantly more time to reach end state (see Figure 48). Figure 48 
highlights one of the challenges of feedback control of flexible systems: it makes 
a difference where measurements are taken. 
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Figure 48. Computed torque control of a flexible DGM using control law 
(5.18): (a) angular displacement; (b) angular rates; (c) control 

torques. 


C. OPTIMIZING THE MANEUVER 


Because the flexible DGM system contains nonlinear terms, the state- 
space representation of the equations of motion takes on a different form. 
Nonetheless, these new dynamics can be used to build an optimal control 
problem to reduce the slew time. From (5.13) and (5.14), the state vector can be 


assembled asy = 


0, 


0 2 co l 


co 2 


01 02 01 02 


T 

, with rotor control vector 


u = 


Ti 


*2 


T 

. For a minimum-time slew, the resulting optimal control problem 


is 
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Minimize 
Subject to 


y i = 3^3 

y 2 = y 4 

• _ 4y 3 y 4 sin v 2 cosy, - y 3 + y 7 - y, + y 5 
6 + 2 sin 2 y 2 + 4 cos 2 y 2 
... _ ~ 2 y\ sin y 2 cos y 2 - y 4 + y 8 - y 2 + y 6 
y *~ 2 
% = y 7 


y 6 = y 8 

y 7 = T| + v3-y 7 +yi-y5 
y 8 = * 2 +y 4 -y 8 +y 2 -y6 
(y, (t 0 ) ,y 2 (*o ) ,y 3 (t 0 ) ,y 4 (t 0 )) = (40°,20°,0,0) 

(y 5 (* 0 )>y 6 (* 0 ) >y 7 (t 0 ) ,y 8 (t 0 )) = (40°,20°,o,o) 

(y, (*/)>y 2 (^).y 3 (^)>y 4 ( ? /)) = (o,4o°,o,o) 

(y 5 (^/).y 6 (*/).y 7 ( ? /)^y 8 ( ? /)) = (o,4o°,o,o) 

|yjsi*%,i = 3,4 
\u\ <1.5 N-m, i = 1,2 


(5.19) 


The challenge in working with these dynamics, given in physical space, is 
that the motion is highly coupled, so it is difficult to constrain the flexible effects. 
To do this correctly, the dynamics should be given in modal space. However, we 
cannot develop the modal model of this system in its current form due to the 
presence of the nonlinear terms. Through arithmetic manipulation, it is possible, 
however, to treat the nonlinear terms as a fictitious torque input, 


T = 

T 


4y 3 y 4 siny 2 cosy 2 - 2z 5 sin 2 y 2 - 4z 5 cos 2 y 2 


t 2 


-2y 2 siny 2 cos y. 


This step uncouples the motion of 0, and 0 2 and provides a linear system of the 
familiar form 


M,0 + C(0-0) + K(0-0) = t 

M 2 0-C(0-0)-K(0-0) = t 
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(5.21) 



Using (5.21), modal analysis can now be easily performed. 

Modal analysis of the unforced system produces the following symmetric 
FRF matrix: 


H(S) = 


5 2 + 5 + l 

6s 4 + 7s 3 + 7s 2 
0 

5 + 1 

65 4 + 75 3 + 75 2 
0 


0 

5 2 + 5 + 1 
25 4 + 35 3 + 35 2 

0 

5 + 1 

25 4 + 35 3 + 35 2 


5 + 1 

65 4 + 75 3 + 75 2 
0 

65 2 + 5 + 1 
65 4 + 75 3 + 75 2 

0 


0 

5 + 1 

25 4 + 35 3 + 35 2 
0 

25 +5 + 1 
25 4 + 35 3 + 35 2 


(5.22) 


Referring to (5.22), the system has four degrees of freedom, yet each 
column of the FRF matrix includes only two poles. The FRF matrix shows 
that torque applied to one axis has no effect on the other (recall that 
X 2 = H 2l F x + H 22 F 2 + H 23 F 3 + H 2A F A ). This makes sense since in the modal model, 
the motion of each axis depends only on the torques applied to that axis. Indeed, 
the coupling between the two axes results from the nonlinear terms, which have 
been relabeled as fictitious torques. 

By removing the nonlinear terms it is possible to analyze the motion of two 
independent systems assuming the simple case where f, = f 2 = 0 (the next 
section discusses how to recover the full nonlinear model). The first system is 
characterized by the azimuth terms: 


6 

0 

o, 

+ 

-0i 

+ 

i 

_cd 

i 

_i 


0 

0 

1 

0i 


01-0, 


0i -0i 


T 


The second system is characterized by the elevation terms: 


2 

0 

1 

CN 

:<£> : 

l_ 

+ 

Q 2 < f > 2 

+ 

d 2 02 


0 

0 

1 

02 


02 _ ^2 


1 

CN 

1 

_1 


^2 


(5.23) 


(5.24) 
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In this context, an optimal control problem definition in physical space may thus 
be defined as follows: 


Yi = 

% = 

... _ --^3 + 3^7 - 3^1 + v 5 
^ 5_ 6 

. -> , 4+.y»-.y2+>'t 

^ 4_ 2 

y 5 = yi 

% = y% 

yi = T\+yi-yi + y\-ys ■ (5.25) 

% =^2+^4-^8+^2- ft 

(y, fa) ,y 2 (to )o5 (to ) >y 4 (t Q )) = (40°,20°,o,o) 

(y 5 (t 0 )^6 (to )^7 (to) ^8 (to )) = (40°,20°,0,0) 

(^1 ( ? / )^2( ? /)^3( ? /)^4 (^/)) = (0, 4 0°,0,0) 

(y 5 (*/)>ft (f/)«ft (*/).y 8 ( ? /)) = (o,40°, o,o) 

|y,.|<irad/,/ = 3 ,4 
\uj\ <1.5 N-m, i = 1,2 

Figure 49 shows the results of the maneuver optimization. However, because the 
nonlinear terms are negated, this result is not valid for the nonlinear system. 
Figure 50 illustrates the maneuver that results from implementing the control 
history from (5.25) to control the nonlinear dynamics of the original problem 
definition given by (5.19). Clearly, the nonlinear system does not reach its 
desired endpoint. 


Minimize 
Subject to 
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Figure 49. Optimized maneuver for flexible, double-axis gimbal with 
nonlinear terms removed (c = k = 1): (a) link displacement; (b) 
link angular rates; (c) control torque. 
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Figure 50. Optimal control solution from Figure 49(c) applied to control 
the nonlinear dynamics of (5.19): (a) link displacement; (b) link 

angular rates. 


It is now illustrated how to recover the original nonlinear dynamics while 
operating in the linear modal space. As was shown previously, the nonlinear 
terms of the DGM may be considered as a fictitious torque input by moving them 
to the right side of the dynamics equations. Modal analysis for this revised 
system may be carried out in much the same manner as in previous examples. 
By treating the nonlinear terms as torque inputs, the modal parameters are 
effectively held constant throughout the maneuver even as the antenna’s 
configuration modifies the inertia of the system over time. Instead of having to 
update the modal model at each computational node, the dynamics for optimal 
control are based on a single nominal modal model. The mass matrix 
M = diag [ 6 2 11 ], which produces a constant eigenvector matrix 
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(5.26) 


10 10 
0 10 1 
10-60 
0 10-2 


Matrix (5.26) is used to transform physical to modal coordinates (and vice-versa). 
The associated modal matrices are 


M = diag[ 7 3 42 6 | 

C = K = diag[ 0 0 49 9 j 


(5.27) 


Recall that for proportionally-damped systems the eigenvectors are the same as 
those of the equivalent undamped system; the modal matrices are scaled the 
same as the eigenvectors. These matrices yield the following natural frequencies 
and critical damping ratios: 


co \ 
0 0 3 

co 4 


= co 2 = 0 
V42 


1.0801 


= — « 1.2247 
2 


C, = C 2 = o 
r _ V42 , 
12 

~ V6 


= 0.5401 

0.6124 


(5.28) 


Note that (5.28) provides two rigid modes (the motion of the DGM links) and two 
flexible modes representing the joint flex. The equations of motion in modal 
space are 

1 q i = u[ 

(5.29) 

42<7, + 49<7, + 49g 3 = u' 3 
6q 4 + 9q 4 + 9 q 4 = u 4 
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where u' = V r u . The first and second equations (for their respective modes) 
result from purely rigid-body motion. 


A new state vector y : 


lT 


4l 4 2 43 44 4l 42 43 44 


is now defined 


with control u = [ f t ] where f = 


T, T 2 


and r = 


T, T 2 


The 


optimal control problem for the nonlinear system expressed in terms of the fixed 
modal model with fictitious torque inputs is given as 


Minimize 
Subject to 


J = t f 

y. = y 5 
y 2 = y 6 
y 3 = y 7 
y 4 = y 8 

y s = y( u i+ u 3 ) 
y6=^( U 2+ U 4) 

117 7 

v 7 = — u, - tin — y 7 —v, 

7 42 1 7 3 6' 7 6' 3 


1 1 


y 8 ^ u 2 


y*-~y* 


6 3 2 2 
(yi( ? o),y2( ? o),y 3 (io),y4( ? o)) = V' 1 [ 

(y 5 ( ? o),y6( ? o),y7(io),y8( f o)) = v _1 [ 

(vi(r / ),y 2 (o),^3(o)’T 4 (? / )) = V“'[ 

( v 5 ) , V 6 ( V ) ,3^7 (0 ) ^8 ( ? / )) = V_1 

U-|<V'[ 1 1 1 1 J ra %,/ = 5,6,7,8 

MO-a (01 


40' 

3 

20° 40° 

20 " J 



i T 


0 

0 

0 0 J 


0 

40° ±180° 

±180° 

0 

0 

+180 y 

/sec 

± 180 / 


sec 


^2 (0 T>2 (l 


\u 3 [t 


\u 4 {t 


< 


0 

0 

1.5 N-m 
1.5 N-m 


(5.30) 
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In solving (5.30), the torque limits are applied in modal coordinates through the 
transformation in the dynamics and the angular rate limit at each node as path 
constraints. Path constraints are also used to ensure that fictitious torques f 
admit the correct nonlinear effect. It is the presence of these two path constraints 
that allows the full nonlinear dynamics to be incorporated with the fixed modal 
model. This aspect represents a new innovation of this thesis. 

As in Chapter IV, we first study the behavior of the system without any 
constraints on flexibility. Figure 51 shows the optimized maneuver resulting from 
this approach after transformation back into physical coordinates. The azimuth 
axis rotates the greater angular distance, requiring the greater control torque, but 
the controller is able to take advantage of the added freedom in the elevation 
axis. The torque profiles of Figure 51(c) are quite different than those 
corresponding to computed torque control. This allows the maneuver time to be 
reduced. 
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Figure 51. Optimal slew for modal system with fictitious torques: (a) axes 
displacement; (b) axes angular rates; (c) input torques. 

Since we are only concerned with the torques supplied by the spacecraft 
actuators, only t 1 and t 2 are needed to implement the maneuver. The values of 
the fictitious controls are, however, included for reference (see Figure 52). 
Applying t, and t 2 to control the full nonlinear dynamic model produces the 
correct maneuver, as shown in the first part of Figure 53 (to the left of the vertical 
dashed line). Thus, the new optimal control problem formulation allows the 
nonlinear effects to be properly recovered. 
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Figure 52. Fictitious control, f. 


In solving the optimal control problem, the vibrations have not been 
constrained in any way. Figure 53 also shows the behavior of the system after 
the end of the minimum-time slew (to the right of the vertical dashed line). 
Because of the energy stored by the “spring” in the joint, each rotor and axis 
continue to spin. This maneuver is not feasible for operations since it is required 
for the antenna to remain pointed at its target. 
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Figure 53. Flexible spillover for minimum-time slew of flexible double 

gimbal: (a) axis; (b) rotor. 


To complete the maneuver with zero flexible motion, the endpoint 
constraints in (5.30) can be modified to 

(-T( ? /)^2(^)^3(^)^4^/)) = V r [ 0 40° 0 40° ] 

(5.31) 

(^ 5 ( f r)) = [ o 0 0 0 ] =( 0 , 0 , 0 , 0 ) 

This ensures that the rotor and link angles match at the end of the maneuver so 
that there is no residual energy stored in the spring. One would expect the slew 
time to increase because the controller must now force the total motion of the 
system to zero. Figure 54 illustrates that this conjecture is correct, while Figure 
55 verifies that motion ceases at t f . Note also that in the zero flex case the 

flexible-body modes q 3 and q 4 are excited (see Figure 56) because they are not 
otherwise constrained. 
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Figure 54. Minimum-time maneuver with endpoint flex constraint, 
t f = 6.496 , k = 1: (a) axes displacement; (b) axes angular rates; 

(c) input torques. 



(a) (b) 

Time (sec) Time (sec) 


Figure 55. Flexible spillover for minimum-time slew of flexible double 
gimbal with terminal flex constraints: (a) axis; (b) rotor. 
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Figure 56. Displacement of system in modal space, zero flex constraints. 


We can easily study the effects of limiting flexible-body motion now that 
we are operating in modal space. In the context of this problem, such a constraint 
would limit the relative displacement between the rotor and the link. To illustrate 
this, the following constraint may be added: 


1^1 < 10 ° 
kl^ 10 ° 


(5.32) 


Figure 57 shows that due to (5.32) the maneuver time increases dramatically 
because the control torque must be modulated to manipulate the relative 
displacement at each joint. 
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Time (sec) 

Figure 57. Minimum-time maneuver with flex constraints (k = 1, t f = 11.51) 

: (a) axes displacement; (b) axes angular rates; (c) input torques. 

Figure 58(a) shows that for this solution the flexible-body modes are 
properly constrained and that the rigid-body modes respond better to the input 
torques since much less energy is stored in the spring of the joints. Figure 58(b) 
confirms that the rotor and axis more closely track each other, showing that 
constraining the flexible modes allows the elastic displacement to be controlled. 
The final displacements of the rotors and links are the same, so the antenna 
does not oscillate at the end of the slew. 
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Figure 58. Flexible effects for maneuver with flex constraints: (a) modal 

space; (b) physical space. 


If the endpoint flex constraint on vibrations were removed, the antenna 
would continue to move after the controller shuts off (see Figure 59). The motion 
is, however, slower than that in Figure 53 because of the imposed constraint on 
the displacements of the flexible modes. Thus, the proper combination of 
constraints on the flexible modes both during the slew maneuver and at the 
terminal time appears to be the best way to reduce slew times for flexible 
systems. 
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Figure 59. Link angle for maneuver with constrained flexibility. 

A comparison of the slew maneuvers for the DGM in this chapter is shown 
in Table 2. The computed torque maneuver (without the nonlinear terms) 
executes nearly five times slower than the optimal control maneuvers. The rigid 
DGM executes the fastest but does not take into account the flexible effects of 
the system. The optimal control solutions for the flexible DGM provide the most 
realistic representation of the slew maneuvers while enabling designers to set 
constraints on flexible motion. Even with the constraints, these slew maneuvers 
execute faster than those using conventional control. 


Table 2. Comparison of slew times (in seconds) for a flexible DGM. 



Computed 
torque control 

Optimal 

control 

Time 

Savings 

Unconstrained 

20 

4.5 

77.5% 

Zero terminal flex 

Cannot execute 

6.5 

67.5% 

Zero terminal flex 
and path 
constraint 

Cannot execute 

11.5 

42.5% 
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D. SUMMARY 

This chapter demonstrated the application of the optimal control concepts 
developed in Chapter IV to a nonlinear dynamical system having characteristics 
of a real spacecraft antenna. The effects of constraining modal displacement 
both during the maneuver and at the terminal time were discussed. The fastest 
slews, as expected, occurred with unconstrained modal displacement. These, 
however, require control effort after the slew to eliminate the spillover motion of 
the antenna. This chapter showed that the most effective method to constrain 
flexibility is to define a path constraint and a terminal constraint on flexible 
motion. This ensures a time-optimal solution can be produced that limits 
vibrations and ensures the antenna will remain pointed at the desired target after 
the slew is complete. 

The key innovation in Chapter V and this thesis is the separation of the 
nonlinear terms from the model by representing them as a fictitious torque. This 
allowed a nominal modal model to be utilized for optimal control while retaining 
the nonlinear effects as part of the optimal control solution. Using this approach, 
the modal properties of the system become readily apparent, and the optimal 
control solution operates as desired when applied to the full nonlinear dynamical 
system in physical space. 
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VI. CONCLUSION 


This thesis addressed a problem in satellite design and operations where 
a control system must accommodate the effects of slewing a flexible appendage 
while maintaining accurate orientation and minimizing disturbances on the 
spacecraft bus and other payloads. The TDRS, for example, has limited 
availability of its single access antennas due to the capabilities of its current 
control systems. Typical slew maneuvers are executed slowly to reduce 
disturbances and avoid exciting flexible modes, thereby reducing customer 
access times. This thesis resolved this issue by determining optimal control effort 
for rapid slew maneuvers by taking advantage a satellite structure’s modal 
properties rather than discounting them. This thesis presented a new way of 
finding optimal controls for flexible systems with complicated nonlinear dynamics. 
Typically, the nonlinearities are simplified or assumed away to make analysis 
easier. However, this often yields inaccurate results, especially for large slew 
maneuvers. 

A simple linear mass-spring-damper system was used to illustrate the 
value of modal analysis to understanding the behavior of flexible structures. 
Modal analysis of a notional spacecraft antenna was then explored by using a 
3-DOF mass-spring-damper system model where applied forces cause coupled 
motion (translation) of each DOF. Computed torque control was used to 
demonstrate that flexible effects can only be constrained by increasing slew time, 
which counters the goal for control system design. 

Pontryagin’s Principle was applied next to the same 3-DOF system, and 
an optimal control problem was defined. Through modal analysis of the system, 
the motion of each DOF was uncoupled in modal space. Flexible effects could 
then be constrained in modal space, so that the optimal control solution could be 
used to reduce slew times while simultaneously avoiding excitation of the modes. 
These results showed that it is possible to apply motive forces in a manner that 

utilizes flexible properties advantageously. 
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Optimal control and modal analysis were applied to a simplified model of 
an actual spacecraft antenna mounted on a DGM. The DGM was updated with 
elastic joints to model flexibility. The nonlinear terms in the equations of motion 
presented a problem for modal analysis in that the modal model was time- 
dependent. A novel approach was used to recast the nonlinearities as fictitious 
torque inputs to the system. This approach removed the time-dependency, yet 
provided a nominal modal model that accurately described the motion of the 
original nonlinear system. It also enabled the application of effective flex 
constraints for optimal control solutions. In the optimal control solution, the 
controller used the available torque to execute a rapid slew while meeting the 
desired flex constraints. The maneuver time for the flexible DGM was decreased 
by approximately 42% compared to the conventional computed torque control 
solution. Thus, the throughput of an antenna system such as the TDRS single 
access antenna could theoretically be improved using the results of this thesis. 

A. TRANSITION TO PRACTICE 

Transforming the dynamics of a system from physical space to modal 
space has other advantages. Realistically, a control system designer would be 
given the natural frequencies and damping ratios of a structure rather than 
having to derive them. If this is the case, optimizing the control is simple. The 
natural frequencies and damping ratios for each mode easily translate to the 
coefficients of the modal space dynamics equations, as demonstrated in 
Chapter II. These coefficients comprise the modal mass, damping, and stiffness 
matrices and are used to determine the modal behavior of the system. This 
model can then be incorporated with a nonlinear dynamics model as was done in 
Chapter V to develop rapid slew maneuvers that accommodate flexible effects 
according to the specification of the mission. 

B. FUTURE WORK 

Future research can explore the modeling of the behavior of systems with 
hysteretic and non-proportional damping. Uncoupling such systems is far more 
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complicated than those presented in this thesis. The simplified DGM used in this 
thesis assumed zero distance between the azimuth and elevation rotors and the 
spacecraft body; future models can incorporate these offsets to study the impact 
of more complicated dynamics on the controller’s ability to minimize vibrations 
and maintain accurate pointing. This model also assumed unit values for stiffness 
and damping for both axes and did not address the impact of different constants 
for each axis. Transition to practice would require working with the correct model 
for the flexible appendage in question. 

This thesis applied elasticity to the gimbal joints of a spacecraft antenna to 
model the flexibility of the system. Another method would be to model the 
behavior of a flexible cantilever attached rigidly to a base structure. The results of 
such a model could then be added to those of a model with elastic joints to 
extend the applicability of the results. 
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